How to solve a system of second order nonlinear differential equations with boundary conditions
4 次查看(过去 30 天)
显示 更早的评论
Hello, I am having troubles solving a system of second order nonlinear equations with boundary conditions using MATALB Here is the equations:
f''(t)=3*f(t)*g(t) ; g''(t)=4f(t)*g(t);
the boundary conditions are: f(0)=1.5 et g'(o)=0; g(2)=3 et f'(2)=f(2)
14 个评论
Paul
2021-6-13
编辑:Paul
2021-6-13
Can someone explain how specifying initial conditions at t = 0 in that call to ode45 will result in the solution satisfying the desired conditions of the solution at t =2? I tried calling ode45 with a tspan = [ 0 2] to see what the solution looks like at t = 2, but ode45 failed.
Also, it looks like the ordering in F and therefore in odeFun is: g(t), g'(t), f(t), f'(t), which is inconsistent with the initial condition vector input to ode45.
Is there a reason why the initial problem statement is in terms of functions f(t) and g(t), but the code uses y, z, x, and theta?
Why are the constants 5 and 7 in the code? They are not in the differential equations in the original question. I thought they were when I first saw the question, but they are not there now.
Alex Sha
2021-6-18
Refer to the result below please:
g(0) = -1.44989398376437
f'(0) = -2.78301182421201
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/657825/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/657830/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/657835/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/657840/image.jpeg)
采纳的回答
J. Alex Lee
2021-6-12
This is a boundary value problem, so probably the most straightforward way is to use bvp4c or bvp5c if you want to solve numerically, and if so don't bother with syms at all. It looks like you already figured how to break your 2nd order equations into first order ones so in addition to f and g you must have j = f' and k = g' or something. So your original question about the BCs you can pose as
0=k(0) and j(2)=f(2)
4 个评论
Paul
2021-6-13
编辑:Paul
2021-6-13
Because of your original answer I looked at bvp4c for the first time and I didn't understand how to set the boundary conditions. Now I do. Thanks for posting this.
Why does your iguess function look the way it does? I mean, it wasn't obvious to me from looking at the equations what the guessed solution should be.
Did you check how closely the solution satisfies the differential equations? I'm not doubting the result, just curious as to how accurate it is.
J. Alex Lee
2021-6-14
believe me, i had some trouble understanding the syntax for bvp4c/5c as well!
It is not at all intuitive to me either what the initial guess should be. bvp4c failed with the trivial initial guess (zeros), so i just picked something random.
I did at least check the BCs are satisfied, but no I didn't really check the ODE...I figured that can be tested by running your IVP integration.
更多回答(2 个)
Paul
2021-6-13
I was not able to solve the problem for the boundary conditions specified at t=2. But I was able to get a resonable looking solution for the same problem with boundry conditions specified at t = 0.5 by posing it as an optimization problem to solve for the unknown initial conditions that result in the boundary conditions being satisfied.
I'm assuming the equations to be solved are:
f''(t) = 3*f(t)*g(t) + 5
g''(t) = 4*g(t)*f(t) + 7
with initial conditions:
f(0) = 1.5, g'(0) = 0
and boundary constraints defined at tf = 0.5:
g(tf) = 3, f(tf) = f'(tf)
Here is the code:
odeFun = @(t,y) ([y(2); 3*y(1)*y(3) + 5; y(4); 4*y(1)*y(3) + 7]);
tspan = [0 .5];
% solve for the unknown initial conditions with an initial guess
x0 = fminsearch(@myfun,[0;0.01]);
% compute the solution
y0 = [1.5 x0(1) x0(2) 0];
[t,y] = ode45(odeFun,tspan,y0,odeset('MaxStep',0.001));
figure;
plot(t,y),legend('f','fprime','g','gprime');
% verify the solution satisfies the boundary conditions
[1.5-y(1,1) 0-y(1,4) y(end,1)-y(end,2)]
% (approximately) verify the solution satisfies the differential equation
for ii = 1:4
ydot(:,ii) = gradient(y(:,ii),t);
yddot(:,ii) = gradient(ydot(:,ii),t);
end
figure
subplot(211);
plot(t,yddot(:,1) - (3*y(:,1).*y(:,3) + 5)),grid
subplot(212);
plot(t,yddot(:,3) - (4*y(:,1).*y(:,3) + 7)),grid
function f = myfun(x)
% y1(t) = f(t)
% y2(t) = f'(t)
% y3(t) = g(t)
% y4(t) = g'(t)
% y1p(t) = f'(t) = y2(t);
% y2p(t) = f''(t) = 3*f(t)*g(t) + 5 = 3*y1(t)*y3(t) + 5
% y3p(t) = g'(t) = y4(t)
% y4p(t) = g''(t) = 4*g(t)*f(t) + 7 = 4*y1(t)*y3(t) + 7
% y0 = [f(0) f'(0) g(0) g'(0)] = [1.5 x(1) x(2) 0]
odeFun = @(t,y) ([y(2); 3*y(1)*y(3) + 5; y(4); 4*y(1)*y(3) + 7]);
tspan = [0 .5];
y0 = [1.5 x(1) x(2) 0];
[t,y] = ode45(odeFun,tspan,y0,odeset('MaxStep',0.001));
f = (y(end,3) - 3)^2 + (y(end,2) - y(end,1))^2;
end
My quick experiments showed the solution blows up pretty quickly for tf > 1, which is probably why ode45 had trouble. I didn't experiment with any of the ode45 parameters or with the initial guess for the initial conditions to see if a solution could be obtained for tf = 2. Maybe that's doable.
12 个评论
Alex Sha
2021-6-20
Hi, Lewis, other than Matlab, the results given above are obtained by a software package 1stOpt, the code is:
Constant q=5;
ODEStep = 0.05;
SubjectTo f'[1]=q*f[1];
Variable t,f,g,g';
ODEFunction
f'' = 3*f*g + 5;
g'' = 4*g*f + 7;
Data;
0,1.5,NAN,0
1,NAN,3,NAN
Running above code in 1stOpt, will produce outcomes like below:
Root of Mean Square Error (RMSE): 2.35047536989441E-11
Sum of Squared Residual: 2.20989378579211E-21
Parameter Best Estimate
-------------------- -------------
g Initial Value -0.326532260011723
f' Initial Value -3.52601348184548
SubjectTo:
f'[1]-(5*f[1]): -2.22044604925031E-16
====== Output Results ======
File: Data File-1
No t Target f Calculated f Target g Calculated g Target g' Calculated g' Target f' Calculated f'
1 0.05 NAN 1.32818918148094 NAN -0.320129119247432 NAN 0.258282359070272 NAN -3.34480171254278
2 0.1 NAN 1.16569586434799 NAN -0.30046930996834 NAN 0.530610920561326 NAN -3.15305529142449
3 0.15 NAN 1.01307875342995 NAN -0.266807892402695 NAN 0.818513829052272 NAN -2.94962811005628
4 0.2 NAN 0.870934501640508 NAN -0.218349329332252 NAN 1.12250210907913 NAN -2.73413690003614
5 0.25 NAN 0.739861393725033 NAN -0.154295907763186 NAN 1.44217044756334 NAN -2.50688564617297
6 0.3 NAN 0.620427266113702 NAN -0.073890512455263 NAN 1.77632399265239 NAN -2.26877048735619
7 0.35 NAN 0.513142785998624 NAN 0.0235477461809992 NAN 2.12313814093348 NAN -2.02115987614537
8 0.4 NAN 0.41844138318349 NAN 0.138596774550519 NAN 2.48035356307457 NAN -1.76574830953955
9 0.45 NAN 0.336667130462941 NAN 0.271715336379485 NAN 2.84550444881187 NAN -1.50438514523658
10 0.5 NAN 0.268071726398618 NAN 0.423239029750086 NAN 3.21617462381618 NAN -1.23888251398334
11 0.55 NAN 0.212821493462363 NAN 0.593389617958113 NAN 3.59027432197821 NAN -0.970807740361826
12 0.6 NAN 0.17101502805802 NAN 0.782298562875354 NAN 3.96633043332803 NAN -0.701265656849458
13 0.65 NAN 0.142711904692132 NAN 0.990045297177201 NAN 4.34378530220257 NAN -0.43067450519355
14 0.7 NAN 0.127972731300827 NAN 1.21671063144516 NAN 4.72330379377781 NAN -0.158535636512123
15 0.75 NAN 0.126910967897541 NAN 1.46244584569714 NAN 5.10709547531324 NAN 0.116808124639453
16 0.8 NAN 0.139757352023235 NAN 1.72755858998777 NAN 5.49926853604065 NAN 0.398437920185011
17 0.85 NAN 0.166938628631905 NAN 2.01261785758903 NAN 5.90624496153632 NAN 0.691170239306758
18 0.9 NAN 0.209173693825101 NAN 2.31858217663632 NAN 6.33728358626256 NAN 1.00194920785144
19 0.95 NAN 0.267592426685248 NAN 2.64695805257288 NAN 6.80518120428412 NAN 1.34037242136761
20 1 NAN 0.34388571319857 3 3.00000000004701 NAN 7.32725606378444 NAN 1.71942856599285
Alex Sha
2021-6-20
Hi, the chart will be given automatically in 1stOpt UI, also you can make chart yourself by using the result data given above. 1stOpt is a commercial software, there is no download linker as my know.
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!