ODE45 and dsolve result discrepency

조회 수: 36 (최근 30일)
Ridwan Hossain
Ridwan Hossain 2015년 8월 10일
댓글: Ridwan Hossain 2015년 8월 13일
I'm having a weird problem. I'm trying to solve a 2nd order ode with both ode45 and dsolve. The results are fine as long as I have non-zero initial conditions but they don't match when the equation has zero initial condition. Any idea why this is happening? Also, which one would be the right choice as I am supposed to implement it in a larger piece of code.
Here is my script:
clc
%clear all
m=3.4e6;
k=3.51e10;
c=13.8e6;
f=7.2578;
[t1,x]=ode45(@pend,[0 5],[0 0] );
j=1;
t2=0:0.01:5;
l=length(t2);
disp2=zeros(l,1);
vel2=zeros(l,1);
for t=0:0.01:5
sol_disp2=dsolve('m*D2x+c*Dx+k*x=f','x(0)=0,Dx(0)=0');
sol_vel2=diff(sol_disp2);
disp2(j)=vpa(subs(sol_disp2));
vel2(j)=vpa(subs(sol_vel2));
j=j+1;
end
subplot(2,2,1)
plot(t1,x(:,1))
subplot(2,2,2)
plot(t1,x(:,2))
subplot(2,2,3)
plot(t2,disp2)
subplot(2,2,4)
plot(t2,vel2)
and the ode45 function:
function dxdt = pend(t,x)
m=3.4e6;
k=3.51e10;
c=13.8e6;
f=7.2578;
x1=x(1);
x2=x(2);
% fun=@(x) sin(x)/z2;
dxdt=[x2; (f-c*x2-k*x1)/m];
end
thanks in advance
  댓글 수: 1
Torsten
Torsten 2015년 8월 11일
Just insert your symbolic solution into
m*D2x+c*Dx+k*x=f, x(0)=x'(0)=0
to see whether it's correct or not.
Best wishes
Torsten.

댓글을 달려면 로그인하십시오.

채택된 답변

Nitin Khola
Nitin Khola 2015년 8월 12일
Hey Ridwad,
I understand you are facing discrepancies in solutions from "dsolve" and "ode45" for zero initial conditions. It appears that the system has faster dynamics compared to the default tolerances in "ode45". You can set the absolute and relative tolerances to smaller values using "odeset" as follows:
>> options = odeset('RelTol', 1e-10, 'AbsTol', 1e-12);
>> [t1,x]=ode45(@pend,[0 5],[0 0],options);
Setting these tolerances to appropriate values get the solutions from the two solvers to match as shown below. Hope this helps.
  댓글 수: 1
Ridwan Hossain
Ridwan Hossain 2015년 8월 13일
Hi Nitin, Thank you very much. It solves the problem.
Ridwan

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기

태그

Community Treasure Hunt

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

Start Hunting!

Translated by