How do I use ode45 for solving this problem and using the solutions as parameters for solving the next system?

1 次查看（过去 30 天）

Stefan Henning 2015-1-19

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/170554-how-do-i-use-ode45-for-solving-this-problem-and-using-the-solutions-as-parameters-for-solving-the-ne

评论： Stefan Henning 2015-1-22

Dear Matlab-community, I got a problem solving some differential-equations-systems.

First I have to say, that I am a beginner in Matlab. I got a system consisting of 6 first order differential equations as a boundary value problem, but I can split them into 3 systems with initial conditions, so I should be able to solve it using ode45. I have to solve the first system and then for solving the second system I need one of the solutions from the first system as a parameter for each x to go on. First I have to solve this system:

(I): Q1’=kappa*Q2 (II): Q2’=-kappa*Q1 (III): -(Q2+dI*kappa)/(E*I)

The (initial)-conditions are: Q1(l)=F1; Q2(l)=F2; Kappa(l)= Ml3/(E*Il), using Il = (B*H^3)/12;

My Matlab code for the function is:

function dy = rigid(x,y)
E = 210000;
I = 3*x; %Funktion für I
dI = 3; %Ableitung dI von I nach dem Weg
dy = zeros(3,1);    % a column vector
dy(1) = y(3)*y(2);
dy(2) = -y(3)*y(1);
dy(3) = (-y(2)+dI*y(3))/(E*I);

and my main program:

clc
clear all 
%Eingangsparameter:
E = 210000; %N/mm²
Ml3 = 264; %Nmm
F1 = 10; %N
F2 = -15; %N
l = 25; %mm     // Länge der Kontur
B = 5; %mm
H = 10; %mm
Il = (B*H^3)/12;
Randbedingung3 = Ml3/(E*Il);
options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);
[X,Y] = ode45(@rigid,[0 l],[F1 F2 Randbedingung3],options);
plot(X,Y(:,1),'-',X,Y(:,2),'-.',X,Y(:,3),'.')

I do not get any values for the solution. I think the problem is that my conditions are not like Q1(0)=…, Q2(0)=…, but at the end of the interval l. I have no idea how to solve this. It would be awesome if someone could pleas help me with this problem and with the following problem:

I need to get the values for y(1), y(2) and y(3) for each x between [0:l]. Then I need to use y(3)(x) which is kappa(x) for solving the next equation:

(i): theta’=kappa with the initial condition theta(x=0)=0 which is easy to solve by getting theta=kappa*x. But therefore I do not know how to get the values for kappa from the solution of the first system for each x.

It would be awesome if someone could please help me out here :-)

17 个评论
显示 15更早的评论隐藏 15更早的评论

Stefan Henning 2015-1-20

在 MATLAB Online 中打开

Thank you Torsten,

I tried using bvp4c, but it does not work because I get some Errors. Here is my code:

function:

function dy = deriv(x,y)
E = 210000;
I = 3*x^2; %Funktion für I
dI = 6*x; %Ableitung dI von I nach dem Weg
dy = zeros(6,1);    % a column vector
dy(1) = y(3)*y(2);
dy(2) = -y(3)*y(1); 
dy(3) = -y(2)+dI*y(3)/(E*I); 
dy(4) = y(3);
dy(5) = cos(y(4))-1; 
dy(6) = sin(y(4));

main:

clc
clear all 
%Eingangsparameter:
E = 210000; %N/mm²
Ml3 = 2640; %Nmm
F1 = 100; %N
F2 = -150; %N
l = 25; %mm     // Länge der Kontur
B = 5; %mm
H = 10; %mm
%boundary conditions:
Il = (B*H^3)/12;
Randbedingung3 = Ml3/(E*Il);
res = @(ya,yb) [ya(4) ya(5) ya(6) yb(1)-F1 yb(2)-F2 yb(3)-Randbedingung3];
solinit = bvpinit(0:l,zeros(1,6)); %zeros setzt alle Annahmen für die Startwerte auf 0 //0.1 gibt die Teilung der Konturlänge in 0.1er Schritten an
sol = bvp4c(@deriv,res,solinit);
plot(sol.x,sol.y(6,:)) %Zahl... in y(...,:) gibt den Graph für die Variable ... der Matrix an // Bsp. y(4,:) entspricht theta3 
plot(X,Y(:,1),'-',X,Y(:,2),'-.',X,Y(:,3),'.',X,Y(:,4),'o',X,Y(:,5),'.o',X,Y(:,6),'-o')
sol.parameters

Errors:

Error using bvp4c (line 251)

Unable to solve the collocation equations -- a singular Jacobian encountered.
Error in main (line 19)
sol = bvp4c(@deriv,res,solinit);

Does anyone know how to solve this problem?

Stefan.

Stefan Henning 2015-1-20

在 MATLAB Online 中打开

Thank you, I tried to fix those 3 aspects, but still get the same error-message. Here is how I changed the code:

function dydx = deriv(x,y)
E = 210000;
I = 3*x; %Funktion für I
dI = 3; %Ableitung dI von I nach dem Weg
dydx = zeros(6,1);    % a column vector
dydx(1) = y(3)*y(2);
dydx(2) = -y(3)*y(1); 
dydx(3) = -y(2)+dI*y(3)/(E*I); 
dydx(4) = y(3);
dydx(5) = cos(y(4))-1; 
dydx(6) = sin(y(4));

sixbc:

function res = sixbc(ya,yb)
%Eingangsparameter:
E = 210000; %N/mm²
Ml3 = 2640; %Nmm
F1 = 100; %N
F2 = -150; %N
B = 5; %mm
H = 10; %mm
%boundary conditions:
Il = (B*H^3)/12;
Randbedingung3 = Ml3/(E*Il);
res = [ya(4); ya(5); ya(6); yb(1)-F1; yb(2)-F2; yb(3)-Randbedingung3];

and main:

clc
clear all 
%Eingangsparameter:
E = 210000; %N/mm²
Ml3 = 2640; %Nmm
F1 = 100; %N
F2 = -150; %N
l = 25; %mm     // Länge der Kontur
B = 5; %mm
H = 10; %mm
%boundary conditions:
Il = (B*H^3)/12;
Randbedingung3 = Ml3/(E*Il);
solinit = bvpinit(linspace(0.1,l), [F1 F2 Randbedingung3 0 0 0]);
sol = bvp4c(@deriv,@sixbc,solinit);
plot(X,Y(:,1),'-',X,Y(:,2),'-.',X,Y(:,3),'.',X,Y(:,4),'o',X,Y(:,5),'.o',X,Y(:,6),'-o')
sol.parameters

I tried to use linspace(0.1, l) to avoid starting at x=0. Do you geht my mistakes? I have to say, I am a novice at matlab.

Best wishes Stefan.

Stefan Henning 2015-1-20

在 MATLAB Online 中打开

Thank you Torsten,

I tried to eliminate all the mistakes I made and finally got a solution without errors. This is why I changed the I and dI into what they really are at my problem. I hope it is correct now, but I got one more question: I would like to get a figure that shows y(1),y(2),...,y(6) over x[0 to l=0.025] How do I get this? Best wishes,

Stefan.

Here is the new code:

main:

clc
clear all 
%Eingangsparameter:
E = 210000000000; %N/m²
Ml3 = 2.640; %Nm
F1 = 100; %N
F2 = -150; %N
l = 0.025; %m     // Länge der Kontur
B = 0.005; %m
H = 0.010; %m
hmin = 0.001; %m
R = 0.0045; %m
%boundary conditions:
Il = (B*H^3)/12;
Randbedingung3 = Ml3/(E*Il);
solinit = bvpinit(linspace(0.1,l), [F1 F2 Randbedingung3 0 0 0]);
sol = bvp4c(@deriv,@sixbc,solinit);
%plot(x,y(:,1),'-',x,y(:,2),'-.',x,y(:,3),'.',x,y(:,4),'o',x,y(:,5),'.o',x,y(:,6),'-o')
plot(sol.x,sol.y(1,:),'-') %,sol.x,sol.y(2,:),'-.',sol.x,sol.y(3,:),'.',sol.x,sol.y(4,:),'o',sol.x,sol.y(5,:),sol.x,sol.y(6,:)) %Zahl... in y(...,:) gibt den Graph für die Variable ... der Matrix an // Bsp. y(4,:) entspricht theta3
sol.parameters

deriv:

function dydx = deriv(x,y)
E = 210000000000; %N/m²
Ml3 = 2.640; %Nm
F1 = 100; %N
F2 = -150; %N
l = 0.025; %m     // Länge der Kontur
B = 0.005; %m
H = 0.010; %m
hmin = 0.001; %m
R = 0.0045; %m
I = (B/12)*((hmin-sqrt(R^2-(x-R)^2)+R)*2)^3; %Funktion für I
dI = (2*B*(x-R)*(hmin-sqrt(R^2-(x-R)^2)+R)^2)/(sqrt(R^2-(x-R)^2)); %Ableitung dI von I nach dem Weg
dydx = zeros(6,1);    % a column vector
dydx(1) = y(3)*y(2);
dydx(2) = -y(3)*y(1); 
dydx(3) = (-y(2)+dI*y(3))/(E*I); 
dydx(4) = y(3);
dydx(5) = cos(y(4))-1; 
dydx(6) = sin(y(4));

sixbc:

function res = sixbc(ya,yb)
%Eingangsparameter:
E = 210000000000; %N/m²
Ml3 = 2.640; %Nm
F1 = 100; %N
F2 = -150; %N
l = 0.025; %m     // Länge der Kontur
B = 0.005; %m
H = 0.010; %m
hmin = 0.001; %m
R = 0.0045; %m
%boundary conditions:
Il = (B*H^3)/12;
Randbedingung3 = Ml3/(E*Il);
res = [ya(4); ya(5); ya(6); yb(1)-F1; yb(2)-F2; yb(3)-Randbedingung3];

Torsten 2015-1-22

Does it not work when you set

y=deval(sol,sol.x);

plot(sol.x,y(1,:),sol.x,y(2,:),sol.x,y(3,:),sol.x,y(4,:),sol.x,y(5,:),sol.x,y(6,:));

Best wishes

Torsten.

Stefan Henning 2015-1-22

Thank you! I did not understand the deval function, or how to use it. Now it works with your code and it gives the exact same results like the plot(sol.x,sol.y(1,:),...) So I guess it does not matter to define yb=L

Thanks a lot

Stefan.

请先登录，再进行评论。

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