solving 4th order ode

19 次查看（过去 30 天）

显示更早的评论

Roi Bar-On 2019-10-1

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/483098-solving-4th-order-ode

评论： Elayarani M 2020-9-23

采纳的回答： darova

在 MATLAB Online 中打开

Hello everyone! I would appreciate your help in this section.

I have this 4th order,non-linear ,homogeneous ode that I would like to solve:

(2B+1+y)y''-k(1+y)y''''=0 with: y(0)=y(L)=0;y'(0)=y'(L)=1; where B,k and L are constants.

This is my matlab code so far:

L=1;B=0.25;k=0.1;
dydt = zeros(4,1);
dydt(1)=y(2);
dydt(2)=y(3);
dydt(3)=y(4);
dydt(4)=((2.*B+1+y(1)).*y(3))./(k.*(1+y(1)));
y1(1)=0; y1(end)=0; y2(1)=1; y2(end)=1;yint=[y1(1);y1(end);y2(1);y2(end)];
tspan = [0 1];
[t,y] = ode45(@(t,y) dydt, tspan, yint);
plot(t,y,'-o')

I think that something is not right. I have almost no expericance with solving ode's with MATLAB, especially with this high order so basically I gathered pieces of codes from things I saw online with the hope that it will assemble to a reasonable solution. I'm not sure I'm in a good path at all. I would be grateful for some guidance. Roi

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

采纳的回答

darova 2019-10-1

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/483098-solving-4th-order-ode#answer_394343

在 MATLAB Online 中打开

Use bvp4c()

You already have written your equation as system of first-order equations

function dydt = ode4(t,y)
    L=1;B=0.25;k=0.1;
    dydt = zeros(4,1);
    dydt(1)=y(2);
    dydt(2)=y(3);
    dydt(3)=y(4);
    dydt(4)=(2.*B+1+y(1)).*y(3))./(k.*(1+y(1));
end

Here is function for your boundary conditions

function res = bc4(y0, y1)
    res = [y0(1)-0        % y(0) = 0
           y0(2)-1        % y'(0) = 1
           y1(1)-0        % y(end) = 0
           y1(2)-1];      % y'(end) = 1
end

See bvp4c

21 个评论
显示 19更早的评论隐藏 19更早的评论

Roi Bar-On 2019-10-2

在 MATLAB Online 中打开

I think I'm missing something. This is the hall code:

function dydx = bvpfcn01(y) % equation to solve

dydx = zeros(4,1);

dydx(1)=y(2);

dydx(2)=y(3);

dydx(3)=y(4);

%L=1;

B=0.25;k=0.1;

dydx(4)=((2.*B+1+y(1)).*y(3))./(k.*(1+y(1)));

end

%--------------------------------

function res = bcfcn(y0, y1)

res = [y0(1)-0 % y(0) = 0

y0(2)-1 % y'(0) = 1

y1(1)-0 % y(end) = 0

y1(2)-1]; % y'(end) = 1

end

%--------------------------------

function main

xmesh = linspace(0,1,50);

solinit = bvpinit(xmesh, [0 0 0 0]);

sol = bvp4c(@bvpfcn01, @bcfcn, solinit);

% y1 = sol.y(3,:);

plot(sol.x, sol.y, '-o')

end

%--------------------------------

I'm getting this error:

Attempted to access y(2); index out of bounds because numel(y)=1.

Error in bvpfcn01 (line 4)

dydx(1)=y(2);

From what I understand MATLAB can't find y(2) for instance because ''y'' is not defined as an array, therefore the numel comment. How can I define it? and when I type in the command window bvpfcn01(y) , should I submit a value of y as an initial guess? Because it doesn't work. Shouldn't it be bvpfcn01(xmesh,y)? I'm really clueless here:( .. Thanks for your patience.

Roi Bar-On 2019-11-8

在 MATLAB Online 中打开

Hi again!

Iv'e noticed that the solution that I'm getting can be only constant - because the initial guess is a scalar and not a variable. I wrote a new code but I was forced to write it down as a bvp with y(0)=0 and y(1)=1, and actually I don't know what's the value of y(1) so I marked that as ''x''. Actually I only know (as presented in the previous code) that y'(0)=y'(1)=0. This is my code:

function shooting_method
clc
clear all
x=0.5;
x1=fzero(@solver,x);
end
function F=solver(x)
options=odeset('RelTol',1e-8,'AbsTol',[1e-8,1e-8]);
[t,u]=ode45(@equation,[0 1],[0 x],options);%ode15s,ode23s are options as well
s=length(t);
F=u(s,1)-1;
figure(1)
plot(t,u(:,1))
xlabel('x')
    ylabel('rho')
    title('S.M')
end
function dy=equation(~,y)
     dy=zeros(2,1);
     dy(1)=y(2);
     A=1;E=1;
     dy(2)=1./A.*log10(E.*y(1));
    
end

I've simplified the problem now and got a more basical ode as can be seen. I marked y(1)=x , because I don't this value. I'm getting weird results. I wanted to know how in this code I can incorporate my conditions y'(0)=y'(1)=0 and not the bvp ones. I just don't know how to use bcfcn here.

I would love to get some help.

Thank you,

Roi

Roi Bar-On 2019-11-19

在 MATLAB Online 中打开

Finally I understood your explanation and the rational, sorry for being so ignorent in this topic.

I've manually ran the code step by step and I see that basically nothing really happens. Basically we're giving the initial guess x=1 and the plot is just getting values of 1 eveywhere because 1 is really a solution. It seems that the code doesn't look for any other solutions and therefore I'm not getting any physical results. This is the final code for now:

function shooting_method
clc
clear all
x=1;
x1=fzero(@solver,x);
end
function F=solver(x)
options=odeset('RelTol',1e-8,'AbsTol',[1e-8,1e-8]);
[t,u]=ode45(@equation,[0 1],[x 1e-8],options);%1e-8 is for the function and x is for the derivative%ode15s,ode23s are options as well
s=length(t);
F=u(s,2)-0; %y'(1)=0
figure(1)
plot(t,u(:,1))
xlabel('x')
    ylabel('rho')
    title('S.M')
end
function dy=equation(~,y)
     dy=zeros(2,1);
     dy(1)=y(2);
     A=1;E=1;
     dy(2)=1./A.*log10(E.*y(1));
    
end

I would be glad for some help. Maybe this method isn't stable enough or the problem is stiff? If so do you happen to know any other method that might help me solve that problem?

Thanks,

Roi