ode45 does not solve on the specified time interval. How do I fix this?
显示 更早的评论
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 个评论
Sam
2018-5-28
Torsten
2018-5-29
I get this message. So you should again check your equations and take a look at the results up to the time when the solver gives up.
Sam
2018-5-29
Jan
2018-5-29
The posted equation contains meaningless information. Remember that
12261830142026331516929221358481899097130334070767616
is rounded to a double internally, so it is the same as:
1.22618301420263e+52
If you are working numerically with the standard double precision, the equation can be simplified remarkably. Nevertheless, this is not the source of the problem.
Sam
2018-5-29
Jan
2018-5-29
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
2018-5-30
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.
回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
2018-5-30
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
