ode45 does not solve on the specified time interval. How do I fix this?

3 views (last 30 days)
Hi, I have a system of differential equations that I want to solve. However, when trying to solve this DiffEq with ode45(@(t,y)odefun(t,y),tspan,y0), where I have put my odefun in a different file, it does produce results, but on a different time interval than I want. Namely, it solves for t = [0, 0.46]*10^-6, while I specified it should solve for t = [0,1].
My code: tspan = [0,1];
y0 = [0 110 -250 15];
[t,Xsolved] = ode45(@(t,y)odefun(t,y),tspan,y0);
odefun (very lengthy equations): function dXdt = odefun(t,y)
t_f = 30;
eq1 = -(82809797410786635*tan(log((- 74530005132491619276355244274867242431640625*y(4)^2*pi^2 - 8612356148643476394993094470262377676800000000000*y(4)^2*pi^2*(9/10 - (17592186044416*y(2))/5918609519667053)^(23/5) - 248801399849700440447422519369183536317307289600000000*y(4)^2*pi^2*(9/10 - (17592186044416*y(2))/5918609519667053)^(46/5) - 12261830142026331516929221358481899097130334070767616*y(2)^2*y(3)^2*(9/10 - (17592186044416*y(2))/5918609519667053)^(46/5))^(1/2)/(pi*y(4)*8633076226496069765625i + y(4)*pi*(9/10 - (17592186044416*y(2))/5918609519667053)^(23/5)*498799959753106275696640000i - 110733148343331824175415296*y(2)*y(3)*(9/10 - (17592186044416*y(2))/5918609519667053)^(23/5)))*1i))/(281474976710656*y(2)); eq2 = - (10019119168824577875*y(2)^2*(1361735765217477681/(28823037615171174400000*(9/10 - (17592186044416*y(2))/5918609519667053)^(27/10)) + (8450000000000*(8500259669165361/(9444732965739290427392*(9/10 - (17592186044416*y(2))/5918609519667053)^(23/5)) + 13/250))/(250119*y(2)^4*cos(log((- 74530005132491619276355244274867242431640625*y(4)^2*pi^2 - 8612356148643476394993094470262377676800000000000*y(4)^2*pi^2*(9/10 - (17592186044416*y(2))/5918609519667053)^(23/5) - 248801399849700440447422519369183536317307289600000000*y(4)^2*pi^2*(9/10 - (17592186044416*y(2))/5918609519667053)^(46/5) - 12261830142026331516929221358481899097130334070767616*y(2)^2*y(3)^2*(9/10 - (17592186044416*y(2))/5918609519667053)^(46/5))^(1/2)/(y(4)*pi*8633076226496069765625i + y(4)*pi*(9/10 - (17592186044416*y(2))/5918609519667053)^(23/5)*498799959753106275696640000i - 110733148343331824175415296*y(2)*y(3)*(9/10 - (17592186044416*y(2))/5918609519667053)^(23/5)))*1i)^2) + 361283/50000000))/147573952589676412928 - (133588255584327705*((4332790137498830962146934784000*y(2)^2)/3184539876935769439309088518619 - (3982870920455782400*y(2))/5918609519667053 + 288000)*((1266637395197952*y(2))/29593047598335265 - (19652298123655411864023597056*y(2)^2)/175149693231467319161999868524045 + (397449804563656125325221541480305270980608*y(2)^3)/1036642641726526473848753494572130715682924789385 + (5520653160719109*y(4)*((4332790137498830962146934784000*y(2)^2)/3184539876935769439309088518619 - (3982870920455782400*y(2))/5918609519667053 + 288000))/731834939447705600000 + 33/2))/(590295810358705651712*((574685827824708321884380135181365943857027170818850816*y(2)^4)/557771182528674706092576514838123659160733159500324835031693855 - (1050791949051857975174900787749300236976128*y(2)^3)/1036642641726526473848753494572130715682924789385 + (119770698800860541596490268672*y(2)^2)/175149693231467319161999868524045 - (6350779162034176*y(2))/29593047598335265 + 229/5)); eq3 = 0; eq4 = (30*((1192349413690968375975664624440915812941824*y(2)^2)/1036642641726526473848753494572130715682924789385 - (39304596247310823728047194112*y(2))/175149693231467319161999868524045 + 1266637395197952/29593047598335265)*((1266637395197952*y(2))/29593047598335265 - (19652298123655411864023597056*y(2)^2)/175149693231467319161999868524045 + (397449804563656125325221541480305270980608*y(2)^3)/1036642641726526473848753494572130715682924789385 + (5520653160719109*y(4)*((4332790137498830962146934784000*y(2)^2)/3184539876935769439309088518619 - (3982870920455782400*y(2))/5918609519667053 + 288000))/731834939447705600000 + 33/2))/((574685827824708321884380135181365943857027170818850816*y(2)^4)/557771182528674706092576514838123659160733159500324835031693855 - (1050791949051857975174900787749300236976128*y(2)^3)/1036642641726526473848753494572130715682924789385 + (119770698800860541596490268672*y(2)^2)/175149693231467319161999868524045 - (6350779162034176*y(2))/29593047598335265 + 229/5) + (16561959482157327*y(4)*(300*y(2)^2*(36766865660871897387/(96970498370224996352000000*(9/10 - (17592186044416*y(2))/5918609519667053)^(37/10)) - (33800000000000*(8500259669165361/(9444732965739290427392*(9/10 - (17592186044416*y(2))/5918609519667053)^(23/5)) + 13/250))/(250119*y(2)^5*cos(log((- 74530005132491619276355244274867242431640625*y(4)^2*pi^2 - 8612356148643476394993094470262377676800000000000*y(4)^2*pi^2*(9/10 - (17592186044416*y(2))/5918609519667053)^(23/5) - 248801399849700440447422519369183536317307289600000000*y(4)^2*pi^2*(9/10 - (17592186044416*y(2))/5918609519667053)^(46/5) - 12261830142026331516929221358481899097130334070767616*y(2)^2*y(3)^2*(9/10 - (17592186044416*y(2))/5918609519667053)^(46/5))^(1/2)/(y(4)*pi*8633076226496069765625i + y(4)*pi*(9/10 - (17592186044416*y(2))/5918609519667053)^(23/5)*498799959753106275696640000i - 110733148343331824175415296*y(2)*y(3)*(9/10 - (17592186044416*y(2))/5918609519667053)^(23/5)))*1i)^2) + 35851247107254511943359375/(86237027846621531955789824*y(2)^4*cos(log((- 74530005132491619276355244274867242431640625*y(4)^2*pi^2 - 8612356148643476394993094470262377676800000000000*y(4)^2*pi^2*(9/10 - (17592186044416*y(2))/5918609519667053)^(23/5) - 248801399849700440447422519369183536317307289600000000*y(4)^2*pi^2*(9/10 - (17592186044416*y(2))/5918609519667053)^(46/5) - 12261830142026331516929221358481899097130334070767616*y(2)^2*y(3)^2*(9/10 - (17592186044416*y(2))/5918609519667053)^(46/5))^(1/2)/(y(4)*pi*8633076226496069765625i + y(4)*pi*(9/10 - (17592186044416*y(2))/5918609519667053)^(23/5)*498799959753106275696640000i - 110733148343331824175415296*y(2)*y(3)*(9/10 - (17592186044416*y(2))/5918609519667053)^(23/5)))*1i)^2*(9/10 - (17592186044416*y(2))/5918609519667053)^(28/5))) + 600*y(2)*(1361735765217477681/(28823037615171174400000*(9/10 - (17592186044416*y(2))/5918609519667053)^(27/10)) + (8450000000000*(8500259669165361/(9444732965739290427392*(9/10 - (17592186044416*y(2))/5918609519667053)^(23/5)) + 13/250))/(250119*y(2)^4*cos(log((- 74530005132491619276355244274867242431640625*y(4)^2*pi^2 - 8612356148643476394993094470262377676800000000000*y(4)^2*pi^2*(9/10 - (17592186044416*y(2))/5918609519667053)^(23/5) - 248801399849700440447422519369183536317307289600000000*y(4)^2*pi^2*(9/10 - (17592186044416*y(2))/5918609519667053)^(46/5) - 12261830142026331516929221358481899097130334070767616*y(2)^2*y(3)^2*(9/10 - (17592186044416*y(2))/5918609519667053)^(46/5))^(1/2)/(y(4)*pi*8633076226496069765625i + y(4)*pi*(9/10 - (17592186044416*y(2))/5918609519667053)^(23/5)*498799959753106275696640000i - 110733148343331824175415296*y(2)*y(3)*(9/10 - (17592186044416*y(2))/5918609519667053)^(23/5)))*1i)^2) + 361283/50000000) + (((8665580274997661924293869568000*y(2))/3184539876935769439309088518619 - 3982870920455782400/5918609519667053)*((1266637395197952*y(2))/29593047598335265 - (19652298123655411864023597056*y(2)^2)/175149693231467319161999868524045 + (397449804563656125325221541480305270980608*y(2)^3)/1036642641726526473848753494572130715682924789385 + (5520653160719109*y(4)*((4332790137498830962146934784000*y(2)^2)/3184539876935769439309088518619 - (3982870920455782400*y(2))/5918609519667053 + 288000))/731834939447705600000 + 33/2))/((574685827824708321884380135181365943857027170818850816*y(2)^4)/557771182528674706092576514838123659160733159500324835031693855 - (1050791949051857975174900787749300236976128*y(2)^3)/1036642641726526473848753494572130715682924789385 + (119770698800860541596490268672*y(2)^2)/175149693231467319161999868524045 - (6350779162034176*y(2))/29593047598335265 + 229/5)))/73183493944770560000 - (30*((1149371655649416643768760270362731887714054341637701632*y(2)^3)/557771182528674706092576514838123659160733159500324835031693855 - (1576187923577786962762351181623950355464192*y(2)^2)/1036642641726526473848753494572130715682924789385 + (119770698800860541596490268672*y(2))/175149693231467319161999868524045 - 3175389581017088/29593047598335265)*((1266637395197952*y(2))/29593047598335265 - (19652298123655411864023597056*y(2)^2)/175149693231467319161999868524045 + (397449804563656125325221541480305270980608*y(2)^3)/1036642641726526473848753494572130715682924789385 + (5520653160719109*y(4)*((4332790137498830962146934784000*y(2)^2)/3184539876935769439309088518619 - (3982870920455782400*y(2))/5918609519667053 + 288000))/731834939447705600000 + 33/2)^2)/((574685827824708321884380135181365943857027170818850816*y(2)^4)/557771182528674706092576514838123659160733159500324835031693855 - (1050791949051857975174900787749300236976128*y(2)^3)/1036642641726526473848753494572130715682924789385 + (119770698800860541596490268672*y(2)^2)/175149693231467319161999868524045 - (6350779162034176*y(2))/29593047598335265 + 229/5)^2 - (82809797410786635*y(3)*tan(log((- 74530005132491619276355244274867242431640625*y(4)^2*pi^2 - 8612356148643476394993094470262377676800000000000*y(4)^2*pi^2*(9/10 - (17592186044416*y(2))/5918609519667053)^(23/5) - 248801399849700440447422519369183536317307289600000000*y(4)^2*pi^2*(9/10 - (17592186044416*y(2))/5918609519667053)^(46/5) - 12261830142026331516929221358481899097130334070767616*y(2)^2*y(3)^2*(9/10 - (17592186044416*y(2))/5918609519667053)^(46/5))^(1/2)/(pi*y(4)*8633076226496069765625i + y(4)*pi*(9/10 - (17592186044416*y(2))/5918609519667053)^(23/5)*498799959753106275696640000i - 110733148343331824175415296*y(2)*y(3)*(9/10 - (17592186044416*y(2))/5918609519667053)^(23/5)))*1i))/(281474976710656*y(2)^2) - (96559323064259661442131689472*y(2))/35029938646293463832399973704809 - 1636073302130688/5918609519667053; dXdt = zeros(4,1); dXdt(1) = eq1*t_f; dXdt(2) = eq2*t_f; dXdt(3) = eq3*t_f; dXdt(4) = eq4*t_f; end
Thanks in advance!
  8 Comments
Jan
Jan on 29 May 2018
I suggest to solve the actual problem of the integration interval at first. But then a massive simplification of the equation should be the next step, if runtime matters.
Are Mjaavatten
Are Mjaavatten on 30 May 2018
Your system is numerically highly unstable, and the solution depends heavily on the integration method and accuracy parameters. Using ode15s, which is probably more robust than ode45 in this case, the solution seems to have a near singularity at around t = 4.4e-6, where y(4) gets very close to zero.
Try
[t,u] = ode15s(@deWringer,[0,1e-5],y0);
plot(t,u,'.')
(I saved your odefun as deWringer.m.)
Most likely, there is something wrong in your derivation of odefun.

Sign in to comment.

Answers (0)

Tags

Community Treasure Hunt

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

Start Hunting!