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

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 :-)

[Torsten]

Be careful: Since y3 depends on x, the solution to theta'=kappa is not simply theta=kappa*x.

I suggest you solve the equations simultaneously using bvp4c.

Best wishes

Torsten.

[Stefan Henning]

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.

[Torsten]

1. Don't start at x=0 ; you'll get dy(3)=infinity there.

2. As starting guess, choose y1=F1, y2=F2, y3=Randbedingung3, y4=0, y5=0 and y6=0.

3. Turn your res-vector into a column vector.

Best wishes

Torsten.

[Stefan Henning]

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.

[Torsten]

I guess there is an error in your equations.

Solving only the first three equations, the solutions of all components blow up tremendously. Maybe a problem of units ? Best wishes

Torsten.

[Torsten]

dy(3) and dydx(3) are different in both codes.

I can't decide which formulation is the correct one.

Best wishes

Torsten.

[Stefan Henning]

I am sorry, I did some mistakes.

dydx(3) = (-y(2)+dI*y(3))/(E*I);

would be the correct code. I corrected all units to meter instead of millimeters:

E = 210000000000;
Ml3 = 2.640;
F1 = 100; 
F2 = -150; 
B = 0.005; 
H = 0.010;

[Torsten]

I used

http://www4.ncsu.edu/~mtchu/Teaching/Lectures/MA581/colsys.f

and the BVP code works fine for your problem.

I can't explain why the MATLAB error message appears.

Best wishes

Torsten.

[Stefan Henning]

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]

x=linspace(0.1,l);
y=deval(sol,x);
plot(x,y(1,:),x,y(2,:),x,y(3,:),x,y(4,:),x,y(5,:),x,y(6,:));

Best wishes

Torsten.

[Stefan Henning]

Thanks a lot! Finally I got some solutions. One last question: how is the solver knowing that yb=l? or do I have to define anywhere?

Thank you

Stefan.

[Torsten]

You mean when you call the deval function ?

It's not necessary to know that yb=l here - the deval function just interpolates the solution from the x-array of the solution saved in the structure "sol" to the points you specify in your vector "x" in the call to deval.

Best wishes

Torsten.

[Stefan Henning]

That is what I meant, thank you. Because I have some conditions like y(2)(l)=F2, and so on, I thought yb=l has to be defined. But I defined it in the x-array like you said. Then the values should be correct, I will check it with FEM.

Thanks a lot for helping me

Stefan.

[Torsten]

You can call the deval-function with sol.x (the mesh selected by bvp4c). Then you can be sure that y(2,end) will be the solution y2 at x=l.

Best wishes

Torsten.

[Stefan Henning]

How do you mean calling the deval-function with sol.x? I don't understand how to.

Best wishes Stefan.

[Torsten]

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]

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.