2017-03-16 7 views
0

私は、内部に追加のプロットを追加したいというコードがあります。私は2つの関数の相対誤差をプロットしたいと思います。私は、エラーメッセージを表示せず、3つのケース全てを1つのグラフとしてプロットするように、これを配置するループがどこにあるのか不明です。また、私の相対的なエラーコードにはいくつか間違いがあるかもしれませんが、それは問題を引き起こしている可能性があります。MATLABでプロットする

これは相対エラーコード: Image 1

rel_error = (y_exact1 - Y(:,2)')./y_exact; %relative error 
figure() 
plot(T,rel_error,'r') 

これは、私はこれはk値のいずれかの所望の結果である

function ivp1() 
clear;clc;close all; 
t=linspace(0,2.5); 
K=[.02 .1 1.5]; 

for i=1:3 
k =K(i); 

[T,Y] = ode45(@prblm1_fun,t,0); %Solving for the approximate solution to the IVP 

figure() 
plot(T,Y) 
hold on 

y_exact1 = 1/(k^2+pi^2)*(pi*exp(k*t)-pi*cos(pi*t)-k*sin(pi*t)); 
y_exact2 = 1/(2*k)*(exp(k*(t-1))-1) + pi/(k^2+pi^2)*(exp(k*t) + exp(k* (t-1))); 
y_exact3 = 1/2/k*(exp(k*(t-1))-exp(k*(t-2))) + pi/(k^2+pi^2)*(exp(k*t) + exp(k*(t-1))) + 1/2/(k-1)*(exp(k*(t-2)) - exp(t-2)); 

for i=1:length(t) 

if t(i)<1 
plot(t(i),y_exact1(i),'mo') 
hold on 

elseif t(i)<2 
plot(t(i),y_exact2(i),'mo') 
hold on 

else 
plot (t(i),y_exact3(i),'mo') 
hold on 

end 
end 
end 
function dy= prblm1_fun(t,y) %This is the function of the IVP for  varying values of t 
dy = zeros(1,1); 
if t < 1 
dy(1)= y(1)*k + sin(pi*t); 
elseif t < 2 
dy(1)= y(1)*k + 0.5; 
else 
dy(1)= y(1)*k + exp(t-2)/2; 
end 
end 
end 

にそれを追加する必要がある関数であります

Image 2

+0

あなたはエラーメッセージを取得していますか?希望の結果の図を貼り付けることができますか? –

+0

@ Dev-iLイメージを追加しました。 – Riya

答えて

0

2つの誤った場所:

1)Yは、ode45の結果から100x1の行列です。 Y(:,2)にアクセスしようとすると、境界外のエラーが発生します。

2)y_exactは定義されていません。 y_exact1y_exact2y_exact3しかありません。

サンプルコード:

function ivp1() 
clear;clc;close all; 
t=linspace(0,2.5); 
K=[.02 .1 1.5]; 

for i=1:3 
    k =K(i); 

    [T,Y] = ode45(@prblm1_fun,t,0); %Solving for the approximate solution to the IVP 

    figure() 
    plot(T,Y) 
    hold on 

    y_exact1 = 1/(k^2+pi^2)*(pi*exp(k*t)-pi*cos(pi*t)-k*sin(pi*t)); 
    y_exact2 = 1/(2*k)*(exp(k*(t-1))-1) + pi/(k^2+pi^2)*(exp(k*t) + exp(k* (t-1))); 
    y_exact3 = 1/2/k*(exp(k*(t-1))-exp(k*(t-2))) + pi/(k^2+pi^2)*(exp(k*t) + exp(k*(t-1))) + 1/2/(k-1)*(exp(k*(t-2)) - exp(t-2)); 

    for j=1:length(t) 
     if t(j)<1 
      y_exact = y_exact1(j); 
     elseif t(j)<2 
      y_exact = y_exact2(j); 
     else 
      y_exact = y_exact3(j); 
     end 
     plot(t(j),y_exact,'mo') 
    end 
    rel_error = (y_exact1 - Y(:,1)')./y_exact; %relative error 
    plot(T,rel_error,'r') 
end 
function dy= prblm1_fun(t,y) %This is the function of the IVP for  varying values of t 
dy = zeros(1,1); 
if t < 1 
    dy(1)= y(1)*k + sin(pi*t); 
elseif t < 2 
    dy(1)= y(1)*k + 0.5; 
else 
    dy(1)= y(1)*k + exp(t-2)/2; 
end 
end 
end 

出力例:

enter image description here