MATLAB: What I should do—–second order variable coefficients ODE?

odeode45second order ode

Hi, matlab person,I have some problems about second order variable coefficients ODE.

Thank you very much for any suggestion or help.

equation: x^2 y''+x y'-sin(y)cos(y)+4*k/pi*x sin(y)^2-x^2 sin(y)cos(y)=0
boundary condition: y(0)=0;y'(0)=-0.126

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.

Best Answer

Your initial equation:

x^2 y''+x y'-sin(y)cos(y)+x sin(y)^2-x^2 sin(y)cos(y)=0

does not contain the ‘-4*k*pi’ term that appears in your second one:

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))];

So should

-4*k*pi

-4*k*pi/x^2

instead?

I found no other problems with your conversion of the first equation to your ‘yprime’ matrix expression. You could simplify it a bit by combining the ‘sin(y)*cos(y)’ terms to ‘(1+1/x^2)*sin(y)*cos(y)’ but I doubt that would make much of a performance difference.

I believe ‘ode45’ will avoid the singularity at x=0 on its own, and search for solutions elsewhere. If you want to be certain about this, contact TMW Tech Support and ask them.

Related Solutions

MATLAB: Error using surf X, Y, Z, and C cannot be complex error

this is John BG <mailto:jgb2012@sky.com jgb2012@sky.com>

The problem is not in solving

x = 0:0.8:pi;
A = sin(x).^(1/3)
B = cos(x).^(1/3)

the operator

.^

solves for either real or complex.

The problem lays in the fact that surf doesn't take in complex values.

To solve the long expression of H, since you are dealing with complex numbers, you have to plot 2 surfs, not one: Real and Imaginary, or Modulus and Phase

clc
theta5=[0:0.1:pi];
theta4=[0:0.1:pi];
[X,Y]=meshgrid(theta4,theta5);
H=(((((cos(Y).^2-cos(X).^2-sin(X).^2+cos(X).^2.*sin(Y).^2+sin(X).^2.*sin(Y).^2).*(cos(X).^2.*cos(Y).^2-cos(X).^2*sin(X).^2+cos(Y).^2.*sin(X).^2+cos(X).^2.*cos(Y).^4.*sin(X).^2+2*cos(X).^2.*sin(X).^2.*sin(Y).^2-cos(Y).^2.*sin(X).^2.*sin(Y).^2-cos(X).^2.*sin(X).^2.*sin(Y).^4))./6+(cos(Y).^2-cos(X).^2-sin(X).^2+cos(X).^2.*sin(Y).^2+sin(X).^2.*sin(Y).^2).^3./27+(cos(X).^2.*cos(Y).^2.*sin(X).^2)./2-(cos(X).^2.*cos(Y).^6.*sin(X).^2)./2-(cos(X).^2.*cos(Y).^2.*sin(X).^2.*sin(Y).^2)./2-cos(X).^2.*cos(Y).^4.*sin(X).^2.*sin(Y).^2).^2-((cos(X).^2.*cos(Y).^2)./3-(cos(X).^2.*sin(X).^2)./3+(cos(Y).^2.*sin(X).^2)./3+(cos(Y).^2-cos(X).^2-sin(X).^2+cos(X).^2.*sin(X).^2+sin(X).^2.*sin(Y).^2).^2./9+(cos(X).^2*cos(Y).^4.*sin(X).^2)./3+(2.*cos(X).^2.*sin(X).^2.*sin(Y).^2)./3-(cos(Y).^2.*sin(X).^2.*sin(Y).^2)./3-(cos(X).^2.*sin(X).^2.*sin(Y).^4)./3).^3).^(1/2)-((cos(Y).^2-cos(X).^2-sin(X).^2+cos(X).^2.*sin(Y).^2+sin(X).^2.*sin(Y).^2).*(cos(X).^2.*cos(Y).^2-cos(X).^2.*sin(X).^2+cos(Y).^2.*sin(X).^2+cos(X).^2.*cos(Y).^4.*sin(X).^2+2.*cos(X).^2.*sin(X).^2.*sin(Y).^2-cos(Y).^2.*sin(X).^2.*sin(Y).^2-cos(X).^2.*sin(X).^2.*sin(Y).^4))./6-(cos(Y).^2-cos(X).^2-sin(X).^2+cos(X).^2.*sin(Y).^2+sin(X).^2.*sin(Y).^2).^3./27-(cos(X).^2.*cos(Y).^2.*sin(X).^2)./2+(cos(X).^2.*cos(Y).^6.*sin(X).^2)./2+(cos(X).^2.*cos(Y).^2.*sin(X).^2.*sin(Y).^2)./2+cos(X).^2.*cos(Y).^4.*sin(X).^2.*sin(Y).^2).^(1/3)+cos(X).^2./3-cos(Y).^2./3+sin(X).^2./3+((cos(X).^2.*cos(Y).^2)./3-(cos(X).^2.*sin(X).^2)./3+(cos(Y).^2.*sin(X).^2)./3+(cos(Y).^2-cos(X).^2-sin(X).^2+cos(X).^2.*sin(Y).^2+sin(X).^2.*sin(Y).^2).^2./9+(cos(X).^2.*cos(Y).^4.*sin(X).^2)./3+(2*cos(X).^2.*sin(X).^2.*sin(Y).^2)./3-(cos(Y).^2.*sin(X).^2.*sin(Y).^2)./3-(cos(X).^2.*sin(X).^2.*sin(Y).^4)./3)./(((((cos(Y).^2-cos(X).^2-sin(X).^2+cos(X).^2.*sin(Y).^2+sin(X).^2.*sin(Y).^2).*(cos(X).^2.*cos(Y).^2-cos(X).^2.*sin(X).^2+cos(Y).^2.*sin(X).^2+cos(X).^2.*cos(Y).^4.*sin(X).^2+2.*cos(X).^2.*sin(X).^2.*sin(Y).^2-cos(Y).^2.*sin(X).^2.*sin(X).^2-cos(X).^2.*sin(X).^2.*sin(Y).^4))./6+(cos(Y).^2-cos(X).^2-sin(X).^2+cos(X).^2.*sin(Y).^2+sin(X).^2.*sin(Y).^2).^3./27+(cos(X).^2.*cos(Y).^2.*sin(X).^2)./2-(cos(X).^2.*cos(Y).^6.*sin(X).^2)./2-(cos(X).^2.*cos(Y).^2.*sin(X).^2.*sin(Y).^2)./2-cos(X).^2.*cos(Y).^4.*sin(X).^2.*sin(Y).^2).^2-((cos(X).^2.*cos(Y).^2)./3-(cos(X).^2.*sin(X).^2)./3+(cos(Y).^2.*sin(X).^2)./3+(cos(Y).^2-cos(X).^2-sin(X).^2+cos(X).^2.*sin(Y).^2+sin(X).^2.*sin(Y).^2).^2./9+(cos(X).^2.*cos(Y).^4.*sin(X).^2)./3+(2.*cos(X).^2.*sin(X).^2.*sin(Y).^2)./3-(cos(Y).^2.*sin(X).^2.*sin(Y).^2)./3-(cos(X).^2.*sin(X).^2.*sin(Y).^4)./3).^3).^(1/2)-((cos(Y).^2-cos(X).^2-sin(X).^2+cos(X).^2.*sin(Y).^2+sin(X).^2.*sin(Y).^2).*(cos(X).^2.*cos(Y).^2-cos(X).^2.*sin(X).^2+cos(Y).^2.*sin(X)^2+cos(X).^2.*cos(Y).^4.*sin(X).^2+2.*cos(X).^2.*sin(X).^2.*sin(Y).^2-cos(Y).^2.*sin(X).^2.*sin(Y).^2-cos(X).^2.*sin(X).^2.*sin(Y).^4))./6-(cos(Y).^2-cos(X).^2-sin(X).^2+cos(X).^2.*sin(Y).^2+sin(X).^2.*sin(Y).^2).^3./27-(cos(X).^2.*cos(Y).^2.*sin(X).^2)./2+(cos(X).^2.*cos(Y).^6.*sin(X).^2)./2+(cos(X).^2.*cos(Y).^2.*sin(X).^2.*sin(Y).^2)./2+cos(X).^2.*cos(Y).^4.*sin(X).^2.*sin(Y).^2).^(1/3)-(cos(X).^2.*sin(Y).^2)./3-(sin(X).^2.*sin(Y).^2)./3;
figure(1);surf(X,Y,abs(H))
figure(2);surf(X,Y,angle(H))

if you find this answer useful would you please be so kind to consider marking my answer as Accepted Answer?

To any other reader, if you find this answer useful please consider clicking on the thumbs-up vote link

thanks in advance for time and attention

John BG

<mailto:jgb2012@sky.com jgb2012@sky.com>

MATLAB: Plot function to graph differential equations

Use a different tspan vector.

When I run your code (using random values for ‘A’ and ‘B’, since they are not supplied), I get:

Warning: Failure at t=5.659410e-01.  Unable to meet integration tolerances without reducing the
step size below the smallest value allowed (1.776357e-15) at time t.

With:

tspan = [0 0.5];

I get a plot.

Also, ‘dx(4)’ as posted is missing a couple parentheses. This corrected version works:

dydx(4) = y(3)*sqrt((y(1)*y(1))+ (y(3)*y(3)))-9.81;

Be sure that is the syntax you intended.

Best Answer

Related Solutions

MATLAB: Error using surf X, Y, Z, and C cannot be complex error

MATLAB: Plot function to graph differential equations

Related Question