second order variable coefficients ODE -so confused !!!

2012 7 12
1 个回答
5 次查看（30 天）

0 个投票

Hi matlab person:
Nowdays, I am solving a second order variable coefficients ODE {x^2 y''+x y'-sin(y)cos(y)+x sin(y)^2-x^2 sin(y)cos(y)=0}. I feel there is no analytical solution, so I want to use ODE to get numerical solution.But I don't know how I should solve second order variable coefficients ODE through ODE45.
Any suggestion, any help is ok.
Thanks very much.
equation: x^2 y''+x y'-sin(y)cos(y)+x sin(y)^2-x^2 sin(y)cos(y)=0
boundary condition: y(0)=0;y'(0)=-0.126

4 个评论
显示 2更早的评论 隐藏 2更早的评论

Yash
Yash 2012-7-12
CHECK MY RUNGE KUTTA CODE FOR THAT, MIGHT HELP, ITS FOR TWO VARIABLES
Following matlab functions will be used. function dxdt = f(t,x,y) % main function dxdt =2*y+x+t;
function dydt = g(t,x,y) % main function dydt =sin(x)+exp(t);
function [tout,yout,xout] = rk2variable(f,g,tspan,y0,N,x0) h = (tspan(2)-tspan(1))/N; t = tspan(1); tout = t; y = y0(:); yout = y.'; x = x0(:); xout = x.'; for n=1:N k1 =feval(f,t,x,y); j1=feval(g,t,x,y); k2=feval(f,t+h/2, x + (h/2)*k1, y + (h/2)*j1); j2=feval(g, t + h/2, x + (h/2)*k1, y + (h/2)*j1); k3=feval(f, t + h/2, x + (h/2)*k2, y + (h/2)*j2); j3=feval(g, t + h/2, x + (h/2)*k2, y + (h/2)*j2); k4=feval(f, t + h, x + h* k3, y + h* j3); j4=feval(g, t + h, x + h* k3, y + h* j3); x = x + (h/6)*(k1 + 2*k2 + 2*k3 + k4); y = y + (h/6)*(j1 + 2*j2 + 2*j3 + j4); t = t+h; yout = [yout; y.']; tout = [tout; t]; xout = [xout; x.']; end
This m file will be used to compute the values at different values of h.
clear tmin = 0; tmax = 1; y0 = 0; h =0.01; N = round((tmax-tmin)/h); x0=0; [tout,yout,xout] = rk2variable(@f,@g,[tmin,tmax],y0,N,x0); plot(tout,yout,tout,xout) legend('Y','X')
Yash
Yash 2012-7-12
Following matlab functions will be used.
function dxdt = f(t,x,y) % main function
dxdt =2*y+x+t;
function dydt = g(t,x,y) % main function
dydt =sin(x)+exp(t);
function [tout,yout,xout] = rk2variable(f,g,tspan,y0,N,x0)
h = (tspan(2)-tspan(1))/N;
t = tspan(1); tout = t;
y = y0(:); yout = y.';
x = x0(:); xout = x.';
for n=1:N
k1 =feval(f,t,x,y);
j1=feval(g,t,x,y);
k2=feval(f,t+h/2, x + (h/2)*k1, y + (h/2)*j1);
j2=feval(g, t + h/2, x + (h/2)*k1, y + (h/2)*j1);
k3=feval(f, t + h/2, x + (h/2)*k2, y + (h/2)*j2);
j3=feval(g, t + h/2, x + (h/2)*k2, y + (h/2)*j2);
k4=feval(f, t + h, x + h* k3, y + h* j3);
j4=feval(g, t + h, x + h* k3, y + h* j3);
x = x + (h/6)*(k1 + 2*k2 + 2*k3 + k4);
y = y + (h/6)*(j1 + 2*j2 + 2*j3 + j4);
t = t+h;
yout = [yout; y.']; tout = [tout; t];
xout = [xout; x.'];
end
This m file will be used to compute the values at different values of h.
clear
tmin = 0; tmax = 1; y0 = 0; h =0.01;
N = round((tmax-tmin)/h);
x0=0;
[tout,yout,xout] = rk2variable(@f,@g,[tmin,tmax],y0,N,x0);
plot(tout,yout,tout,xout)
legend('Y','X')
Jan
Jan 2012-7-12
编辑：Jan 2012-7-12
@Yash: I definitely would not hand-code a Runge-Kutta manually, when Matlab offers more powerful ODE functions.
Jan
Jan 2012-7-12
@petter: Please consider, that the term "urgent" is not apropriate, when the answers are craeted by volunteers. Even the "please" does not shift the tone back to a friendly level sufficently.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Jan
Jan 2012-7-12
编辑：Jan 2012-7-12

0 个投票

You can convert the ODE system of order k for a problem of dimension d into a system of order 1 and dimension k*d. Instructions can be found at Wikipedia or in your Numerics script.
Then the documentation of ODE45 will help you to implement this in Matlab, especially the examples on "doc ode45".

1 个评论
显示 -1更早的评论 隐藏 -1更早的评论

petter
petter 2012-7-12
编辑：petter 2012-7-12
Jan simon, thank you very much for your reply and for your suggestion. I have thought this problem several days and asked other persons. This problem comes from the paper which I am preparing.So I used the term"urgent",may be that is not appropriate. sorry for this.
The following is my code, but there are several problem:
first defining function ***********************************************
function yprime = skyrmion(x,y)
k=0.95;
yprime = [y(2); -1/x*y(2)+1/x^2*sin(y(1))*cos(y(1))-4*k/pi* (1/x)*sin(y(1))^2+sin(y(1))*cos(y(1))]; ***********************************************************
Then transfer function
*******************************************************
xspan=[0.01,10];
y0=[-0.126;pi];
[x,y]=ode45('skyrmion',xspan,y0);
plot(x(:),y(:,2))
******************************************************
problem one:
Because The ODE was divided by x^2, so the x span can't use x==0. whether the code can be modified to not divide x^2.
problem two:
the results get by using my code are not consist with some papers.
so I want to know my code is right or wrong.
Thank you again.

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

提问:

petter
petter
2012-7-12

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by