Equation: X”-6x’+13x = t+3sin(t) Initial Value: x(0)=1 t є [0,1] Method: Runge-Kutta II Step Sizes: h=0.1 , h=0.03

1 次查看（过去 30 天）
Tariq Malik
Tariq Malik 2018-3-30
编辑： Walter Roberson 2018-4-8
I want to solve it by the Matlab only. but facing the Problem . Can someone Please help me out?
  3 个评论
显示 1更早的评论
Tariq Malik
Tariq Malik 2018-3-30
编辑：James Tursa 2018-3-30
This is the code that I am using.
intmin=0;
intmax=1;
numnodes=5;
f=@(t,x) t+3*(sin(t));
inival=1;
h=(intmax-intmin)/(numnodes-1);
t=zeros(1,numnodes);
x=zeros(1,numnodes);
X=zeros(1,numnodes);
t(1)=intmin;
x(1)=inival;
X(1)=inival;
for i=2:numnodes;
t(i)=t(i-1)+h;
k1=f(t(i-1),x(i-1));
k2=f(t(i-1)+h/2,x(i-1)+(h/2)*k1);
x(i)= x(i-1)+h*k2;
X(i)=X(i-1)+h*f(t(i-1),X(i-1));
end
figure
plot(t,x,'.-',t,X,'-*')
hold on
syms t x(t) sym(f)
eqn=diff(x,t)==f(t,x(t));
cond1=x(intmin)==inival;
odesol=dsolve(eqn,cond1);
odesolfun=@(t) eval(odesol);
tt=linspace(intmin,intmax,100*(ceil(intmax-intmin)));
xx=odesolfun(tt);
plot(tt,xx,'r')
hold off

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Abraham Boayue
Abraham Boayue 2018-3-31
编辑：Abraham Boayue 2018-3-31
The first thing you need to do is to write the ode as two first order equations and use the code below. You will be required to supply two initial conditions for the 1s order equations. Use the one that you are given plus another of your choice.
function [t,x,y,N] = Runge2_2eqs(f1,f2,to,tfinal,xo,yo,h)
% This function implements the Rk2 method.
t = to;
N = ceil((tfinal-to)/h);
x = zeros(1,N);
y = zeros(1,N) ;
x(1) = xo;
y(1) = yo;
for i = 1:N
t(i+1) = t(i)+h;
Sx1 = f1(t(i),x(i),y(i));
Sy1 = f2(t(i),x(i),y(i));
Sx2 = f1(t(i)+h, x(i)+Sx1*h, y(i)+Sy1*h);
Sy2 = f2(t(i)+h, x(i)+Sx1*h, y(i)+Sy1*h);
x(i+1) = x(i) + h/2*(Sx1+Sx2);
y(i+1) = y(i) + h/2*(Sy1+Sy2);
end
end
This is the mfile.
xo = 1;
yo = 0;
h = [.1 0.03];
to = 0;
tfinal = 20;
M = ceil((tfinal-to)/h(2));
dx1 = @(t,x1,x2) x2;
dx2 = @(t,x1,x2) 6*x2 -13*x1 + t + 3*sin(t);
% When you reduce the equation to two first order, x will be the solution
% of the ode, i.e x'' and y is its derivative, x'.
for i = 1: length(h)
if (i== 1) % This for the case when h = 0.1
[t,x,y,N] = Runge2_2eqs(dx1,dx2,to,tfinal,xo,yo,h(i));
y1 = x;
y2 = y;
else % and for the case when h = 0.03
[t,x,y,N] = Runge2_2eqs(dx1,dx2,to,tfinal,xo,yo,h(1));
x3 = x;
x4 = y;
end
end
t1 = t(1):(t(end)-t(1))/(M-1):t(end);
figure(1);
subplot(121)
plot(t1,y1, '-o')
hold on
plot(t1,y2,'-o')
legend('Dfx1','Dfx2')
title('Solution to two systems of ODEs using RK2, h= 0.1')
xlabel('x')
ylabel('y')
xlim([to tfinal])
grid
subplot(122)
plot(t,x3,'linewidth',2,'color','b')
hold on
plot(t,x4,'linewidth',2,'color','r')
legend('Dfx1','Dfx2')
title('Solution to two systems of ODEs using RK2, h = 0.03')
xlabel('x')
ylabel('y')
xlim([to tfinal])
grid
% Using ode 45 just to prove that the solution with RK2 is correct.
F = @(t,y) [ y(2); (6*y(2) -13*y(1) + t + 3*sin(t)) ];
t0 = 0;
tf = 20;
delta = (tf-t0)/(201-1);
tspan = t0:delta:tf;
ic = [1 0];
[t,y] = ode45(F, tspan, ic);
figure
plot(t,y(:,1),'-o')
hold on
plot(t,y(:,2),'-o')
a = title('Using ode45');
legend('x','x_{prime}');
set(a,'fontsize',14);
a = ylabel('y');
set(a,'Fontsize',14);
a = xlabel('t [0 20]');
set(a,'Fontsize',14);
xlim([t0 tf])
grid
grid minor;
  1 个评论
Tariq Malik
Tariq Malik 2018-4-7
编辑：Walter Roberson 2018-4-8
intmin=0;
intmax=1;
inival1=0;
inival2=0;
numnodes = 10;
t(1) = intmin;
x1(1)= inival1;
x2(1)= inival2;
t =zeros(1,numnodes);
x1=zeros(1,numnodes);
x2=zeros(1,numnodes);
h=(intmax-intmin)/(numnodes-1);
f1 = @(t,x1,x2) x2;
f2 = @(t,x1,x2) 6*x2-13*x1+t+3*sin(t);
numnodes1=300;
a=zeros(1,numnodes1);
b1=zeros(1,numnodes1);
b2=zeros(1,numnodes1);
a(1)=intmin;
b1(1)=inival1;
b2(1)=inival2;
g=(intmax-intmin)/(numnodes1-1);
F1 = @(a,b1,b2) b2;
F2 = @(a,b1,b2) 6*b2-13*b1+a+3*sin(a);
for i = 2:numnodes
t(i) = t(i-1)+h;
k1 = f2(t(i-1),x1(i-1),x2(i-1));
k2 = f2(t(i-1)+h/2,x1(i-1)+(h/2)*k1,x2(i-1)+(h/2)*k1);
k3 = f2(t(i-1)+h/2, x1(i-1)+(h/2)*k2,x2(i-1)+(h/2)*k2);
x1(i) = x1(i-1)+(h/6)*(k1+4*k2+k3);
x2(i) = x2(i-1)+(h/6)*(k1+4*k2+k3);
end
for i = 2:numnodes1
a(i) = a(i-1)+g;
k1 = F2(a(i-1),b1(i-1),b2(i-1));
k2 = F2(a(i-1)+g/2,b1(i-1)+(g/2)*k1,b2(i-1)+(g/2)*k1);
k3 = F2(a(i-1)+g/2, b1(i-1)+(g/2)*k2,b2(i-1)+(g/2)*k2);
b1(i) = b1(i-1)+(g/6)*(k1+4*k2+k3);
b2(i) = b2(i-1)+(g/6)*(k1+4*k2+k3);
end
figure
plot(t,x1,'.-',t,x2,'-*')
hold on
syms t x(t)
sym(x)
eqn=diff(x,t,2) == diff(x,t)-13*x+t+3*sin(t);
cond1==1;
cond2==0;
odesol=dsolve(eqn,cond1,cond2);
odesolfun=@(t) eval(odesol);
tt=linspace(intmin,intmax,100*(ceil(intmax-intmin)));
xx=odesolfun(tt);
legend("Runge-Kutta 2", "0.1", "Exact");
plot(tt,xx,'r')
hold off

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by