How to solve a non-linear equation where each element of the array is a parameterized function?

Dear all,

I need to solve an equation with the following form:

 [f(a,b,c); g(a,b,c); p(a,b,c)] = [0; 0; q(t)]

Ok, it's more like

 [f(a,b,c,da,db,dc,d2a,d2b,d2c); g(a,b,c,da, db,dc,d2a,d2b,d2c); p(a,b,c,da,db,dc,d2a,d2b,d2c);] = [0; 0; q(t)]

(where "da" means "da/dt" and "d2a" means "d2a/dt2", got using fulldiff())

Or, exactly:

A=
 12*d2Gama*cos(Alpha)*cos(Gama)^2 - (9*d2Beta*cos(Gama))/10 - (9*d2Alpha*cos(Beta)*sin(Gama))/10 - 12*d2Gama*cos(Alpha) + 12*d2Gama*cos(Beta)*cos(Gama)^2 - (9*dAlpha^2*cos(Beta)*cos(Gama)*sin(Beta))/10 - 12*dGama^2*cos(Alpha)*cos(Gama)*sin(Gama) - 12*dGama^2*cos(Beta)*cos(Gama)*sin(Gama) + 12*dGama^2*cos(Gama)^2*sin(Alpha)*sin(Beta) + (9*dAlpha*dBeta*sin(Beta)*sin(Gama))/5 + 12*d2Gama*cos(Gama)*sin(Alpha)*sin(Beta)*sin(Gama); (9*d2Alpha*cos(Beta)*cos(Gama))/10 - 12*dGama^2*cos(Beta) - (9*d2Beta*sin(Gama))/10 + 12*d2Gama*sin(Alpha)*sin(Beta) + 12*dGama^2*cos(Alpha)*cos(Gama)^2 + 12*dGama^2*cos(Beta)*cos(Gama)^2 - 12*d2Gama*cos(Gama)^2*sin(Alpha)*sin(Beta) - (9*dAlpha^2*cos(Beta)*sin(Beta)*sin(Gama))/10 - (9*dAlpha*dBeta*cos(Gama)*sin(Beta))/5 + 12*d2Gama*cos(Alpha)*cos(Gama)*sin(Gama) + 12*d2Gama*cos(Beta)*cos(Gama)*sin(Gama) + 12*dGama^2*cos(Gama)*sin(Alpha)*sin(Beta)*sin(Gama); 12*dGama^2*sin(Beta)*sin(Gama) - (9*dBeta^2)/10 - 12*d2Gama*cos(Gama)*sin(Beta) - (9*dAlpha^2*cos(Beta)^2)/10 + 12*dGama^2*cos(Beta)*cos(Gama)*sin(Alpha) + 12*d2Gama*cos(Beta)*sin(Alpha)*sin(Gama)];
B= [0 0 log(t+1)+sin(t)./(t+1)].';
fsolve(A-B)

but let's try to learn the simplest first.

I need to get expressions for a(t), b(t) and c(t).

I believe there's no exact solution, so, it could be an approximated solution, since it returns expressions Alpha(t), Beta(t) and Gama(t).

How do we do Matlab solve it?

Thanks, very much

Answer 1

[Torsten]

Add the equations

dalpha/dt=alphadot
dbeta/dt=betadot
dgamma/dt=gammadot,

to your system, substitute

dalpha by alphadot
dbeta by betadot
dGama by gammadot
d2alpha by dalphadot/dt
d2beta by dbetadot/dt
d2Gama by dgammadot/dt

in your equations and use ODE15i to solve.

Best wishes

Torsten.

[Marlon Saveri Silva]

Thanks very much Torsten.

Using what you suggest, I made the code below to learn about ODE15i considering a f(t,y,y')=0 problem. I'll let it here if One someday has the same question.

However, could you help me how to include an f(t,x,x',x",y,y',y",z,z',z")=0 inside an ODE15i()? Is it possible?

function res = myfunc(t,y,yp)
%myfunc  Function to test ODE
%
%   See also ODE15I.

%   Marlon
%   Got from Matlab Help.

res = t*y^2 * yp^3 - y^3 * yp^2 + t*(t^2 + 1)*yp - t^2 * y;

---

%%Test ODE15i
clear
clc
syms t y
t0 = 1;                                         %Initial guess to t0
y0 = sqrt(3/2);                                 %Initial guess to y0
yp0 = 0;                                        %Initial guess to yp0
[y0,yp0] = decic(@myfunc,t0,y0,1,yp0,0);        %"calculated guess" to y0 and yp0
odeset('refine',[10]);
[t,y] = ode15i(@myfunc,[1 10],y0,yp0);          %Numerical solution y(t)
ytrue = sqrt(t.^2 + 0.5);                       %Real solution (to compare)
n=4;                                            %Polynomial degree
p=polyfit(t,y, n);                              %Polynomial coefficients
f = polyval(p,t);                               %Values of the polynomial

syms tsim ytrue
yt=0;                                           

for i=1:n+1
    yt=yt+p(i)*tsim^(n+1-i);                %Polynomial wrote as symbolic
end

%Prepare to plot
numb=1:1:t(end);
yt_plot=subs(yt,tsim,numb);

ytp=diff(yt,tsim);
ytp_plot=subs(ytp,tsim,numb);

%plot and compare
plot(t,y,'o',t,f,'r', t, ytrue,'k', numb, yt_plot, 'm',numb,ytp_plot, 'cyan'); hold;
legend('NumericSolution','Polynomial (numeric)','ExactSol','Polynomial (function)', 'yp');

Answer 2

[Marlon Saveri Silva]

For example. I tried this, but it's not running:

function res = myfunc2(t,y,yp)
%myfunc2  Function to test ODE with 2 variables and their diffs
%
%   See also ODE15I.

%   Marlon
%   Got from Matlab Help.

res = [t*y(1)^2 * yp(1)^3 - y(1)^3 * yp(1)^2 + t*(t^2 + 1)*yp(1) - t^2 * y(1); t*y(2)^2 * yp(2)^3 - y(2)^3 * yp(2)^2 + t*(t^2 + 1)*yp(2) - t^2 * y(2)];

---

clear
clc
t0 = 1;                                         %Initial guess to t0
y1_0 = sqrt(3/2);                                 %Initial guess to y1_0
yp1_0 = 0;                                        %Initial guess to yp1_0
y2_0 = sqrt(3/2);                                 %Initial guess to y2_0
yp2_0 = 0;                                        %Initial guess to yp2_0
[t,y] = ode45(@myfunc2,[1 10],[y1_0 y2_0],[yp1_0 yp2_0]);          %Numerical solution y(t)

[Torsten]

I already included the way how to deal with 2nd order derivatives in my reply .

If you have the equation

f(t,y,y',y'')=0,

replace it by the two equations

y'-y1=0;
f(t,y,y1,y1')=0

Best wishes

Torsten.