ode45 does not solve on the specified time interval. How do I fix this?
3 views (last 30 days)
Show older comments
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
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
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.
Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!