Il sistema è il seguente
$ w'''+1/rw''-1/r^2w'=12/h^2w'*( u'+vu/r+1/2*(w')^2 )+(qr)/(2*D) $
$ u''+1/ru'-u/r=-((1-v)/(2r))(w')^2-w'w'' $
Dove:
* r è la variabile dipendente
* h=1 [millimetri]
* E=200000 [MegaPascal]
* v=0.3
* $ D=(Eh^3)/(12(1-v^2)) $
* q=20 [MegaPascal]
(Il sistema fisicamente rappresenta la deflessione di una piastra circolare di acciaio di raggio a e spessore h soggetta a un carico uniformemente distribuito q; w e u sarebbero rispettivamente la deflessione e lo spostamento radiale della piastra in funzione della posizione radiale r )
L'integrazione deve avvenire sull'intervallo [a b], dove:
*a=12
*b=0
Con le seguenti condizioni al contorno:
*w(r=a)=0
*w'(r=a)=0
*w'(r=0)=0
*u(r=a)=0
*u'(r=0)=0
(Tali condizioni fisicamente rappresentano un'incastro.)
Purtroppo le condizioni al contorno non sono definite in un unico punto, ma in due. Per risolvere la cosa ho pensato di utilizzare il metodo di shooting trasformando il problema in un sistema di ODE. Ho quindi prima testato con successo un codice di calcolo in linguaggio Matlab per una singola equazione (in particolare ho risolto un equazione di una discussione già riportata nel forum (https://www.matematicamente.it/forum/vi ... p?t=189081) e poi adattato il codice per questo sistema. Il codice funziona e calcola correttamente (ho confrontato con analisi FEM) ma non se h diventa troppo piccolo (per esempio h=0.1) o se q diventa troppo grande perchè la soluzione oscilla. Mentre invece usando la funzione di matlab per i problemi con i valori ai bordi "bvp4c" si ottiene la soluzione.
Volevo implementare un'altro schema (ad esempio uno schema di shooting multiplo) per ovviare a questo inconveniente. Ho provato a implementare uno shooting multiplo ma non ho ben capito come costruirlo. Chiedo quindi gentilmente se qualcuno può aiutarmi. Allego anche gli script che ho steso. Grazie a tutti un saluto
- Codice:
function shooting_method
clc
clear all
close all
global a h E v q D
a=12;
h=0.1;
E=200000;
v=0.3;
q=20;
D=E*h^3/12/(1-v^2);
x=[0 0];
options=optimset('Display','iter');
x1=fsolve(@solver,x)
end
function F=solver(x)
global a h E v q D
%options=odeset('RelTol', 1e-08, 'AbsTol', [1e-08 1e-08 1e-08 1e-08 1e-08]);
[t,u]=ode(@equation, [a 0.01], [0 0 x(1) 0 x(2)]);
s=length(t);
F=[ abs( u(s,2) ) , abs( u(s,4) ) ]
clf
plot(t,u(:,1),'-or')
xlabel(' r[mm] ');
ylabel(' dy [mm] ');
end
function dy=equation(r,y)
global a h E v q D
dy=zeros(5,1);
dy(1)=y(2);
dy(2)=y(3);
dy(3)=-y(3)/r+y(2)/r^2+12/h^2*y(2)*( y(5)+v*y(4)/r+0.5*y(2)^2 )+q*r/2/D;
dy(4)=y(5);
dy(5)=-y(5)/r+y(4)/r^2-(1-v)/2/r*y(2)^2-y(2)*y(3);
end
Codice per (https://www.matematicamente.it/forum/vi ... p?t=189081)
- Codice:
function shooting_method
clc
clear all
x=0.5;
display('La condizione iniziale equivalente vale')
x1=fzero(@solver,x)
end
function F=solver(x)
options=odeset('RelTol', 1e-08, 'AbsTol', [1e-08, 1e-08]);
a=1;
b=3;
y0=[17 x];
[t,u]=ode45(@equation, [a b], y0, options);
s=length(t);
F=u(s,1)-43/3;
figure(1)
clf
plot(t,u(:,1),'-*r') % plot della soluzione numerica
hold on
plot(t,(t.^2+16./t),'-b') % plot della soluzione esatta
end
function dy=equation(t,y)
dy=zeros(2,1);
dy(1)=y(2);
dy(2)=1/8.*( 32+2.*t.^3-y(1).*y(2) );
end