ode45 does not solve on the specified time interval. How do I fix this?
2 ビュー (過去 30 日間)
古いコメントを表示
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 件のコメント
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 件)
参考
カテゴリ
Help Center および File Exchange で Ordinary Differential Equations についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!