obtaining the values for 'a' using newton raphson for different values of g and plotting the graph.

1 次查看（过去 30 天）
SATYAJIT BHATTACHARJEE
编辑： VBBV 2020-9-30
Hi all,
In the below mentioned code, I tried to get different values of 'a' for a range of 'g'.
Here, Initially I have two 4th order ODEs for vz and uz which are solved using ode 45.
Now the obtained values are put in the equations to get the matrix elements of 4X4 matrix.
The condition for getting consistant value of 'a' is that, the determinant of 4X4 matrix should be zero.
So, initially 'a0' is supplied and the value is to iterated using newton raphson method till convergence.
But, this code is showing some error for "INDEX EXCEEDING THE NUMBER OF ARRAY ELEMENTS"..
PLease help.................
function newtonraphson
a0= -0.5 %input('enter initial guess for complex velocity of perturbations(c)')
itermax= 3000 %input('enter max. iteration')
h=0.1 %input('enter value for stepsize(h)')
tolerance= 10^-10 %input('enter the tolerance limit')
glist=linspace(1,10,10); % Enter the values of g
for m=1:length(glist)
g=glist(m)
for j=1:itermax
F=f(a0,g)
Fprime=(f(a0+h,g)-f(a0,g))/h
a(j)=a0-(F/Fprime)
Fsolution=f(a(j),g)
if abs(Fsolution)<=tolerance
a(j)
break
else
a0=a(j)
end
end
alist(m)=a(end)
end
figure(1); plot(glist, real(alist), 'b')
xlabel('gamma')
ylabel('a-real')
title('plot of a-real vs gamma')
figure(2); plot(glist, imag(alist), 'r')
xlabel('gamma')
ylabel('a-imaginay')
title('plot of a-imaginary vs gamma')
end
function F=f(a,g)
syms uz(z) vz(z) Y
duz=diff(uz);
d2uz=diff(uz,2);
d3uz=diff(uz,3);
d4uz=diff(uz,4);
dvz=diff(vz);
d2vz=diff(vz,2);
d3vz=diff(vz,3);
d4vz=diff(vz,4);
er=0;% relative viscosity eta_r
H=0.1;% Ratio of thickness of gel to fluid
k=1;% gemma (dimensionless velocity)
re=0;% reynolds number
eqn1= (-k^2*d2uz + k^4*uz + i*k^3*a*er*d2uz - i*k^5*a*er*uz + g^2*k^4*uz - 2*i*k^3*g*duz - (re/g)*k^4*a^2*uz + d4uz - k^2*d2uz - i*k*a*er*d4uz + i*k^3*a*er*d2uz - g^2*k^2*d2uz + 2*i*k*g*d3uz + (re/g)*k^2*a^2*d2uz) ==0;
[UF] = odeToVectorField(eqn1);
coupled1 = matlabFunction(UF, 'Vars',{z,Y});
output1=coupled1;
uz10=0;
duz10=0;
d2uz10=1;
d3uz10=0;
zspan1=linspace(-H,0,10); % range of z
z10 = [uz10,duz10,d2uz10,d3uz10]; % array of initial values
[r1,v1] = ode45(@(z,Y)coupled1(z,Y),zspan1,z10)
%% Declaring the 2nd set independent solutions for displacement field at z=-10
uz20=0;
duz20=0;
d2uz20=0;
d3uz20=1;
zspan2=linspace(-H,0,10); % range of z
z20 = [uz20,duz20,d2uz20,d3uz20]; % array of initial values
[r2,v2] = ode45(@(z,Y)coupled1(z,Y),zspan2,z20)
eqn2= d4vz -2*k^2*d2vz + k^4*vz -i*re*k*z*d2vz - i*k*re*dvz + i*k*(re/g)*a*d2vz - re*d2vz + i*k^3*re*z*vz - (re/g)*i*k^3*a*vz;
[VF] = odeToVectorField(eqn2);
coupled2 = matlabFunction(VF, 'Vars',{z,Y});
output2=coupled2;
vz10=0;
dvz10=0;
d2vz10=1;
d3vz10=0;
Zspan1=linspace(1,0,10); % range of z
Z10 = [vz10,dvz10,d2vz10,d3vz10]; % array of initial values
[R1,V1] = ode45(@(z,Y)coupled2(z,Y),Zspan1,Z10)
%% Declaring the 2nd set independent solutions for velocity field at z=1
vz20=0;
dvz20=0;
d2vz20=0;
d3vz20=1;
Zspan2=linspace(1,0,10); % range of z
Z20 = [vz20,dvz20,d2vz20,d3vz20]; % array of initial values
[R2,V2] = ode45(@(z,Y)coupled2(z,Y),Zspan2,Z20)
v1_1=v1(:,1)
v1_2=v1(:,2)
v1_3=v1(:,3)
v1_4=v1(:,4)
v2_1=v2(:,1)
v2_2=v2(:,2)
v2_3=v2(:,3)
v2_4=v2(:,4)
V1_1=V1(:,1)
V1_2=V1(:,2)
V1_3=V1(:,3)
V1_4=V1(:,4)
V2_1=V2(:,1)
V2_2=V2(:,2)
V2_3=V2(:,3)
V2_4=V2(:,4)
uz1=v1_1(end)
duz1=v1_2(end)
d2uz1=v1_3(end)
d3uz1=v1_4(end)
uz2=v2_1(end)
duz2=v2_2(end)
d2uz2=v2_3(end)
d3uz2=v2_4(end)
vz1=V2_1(end)
dvz1=V2_2(end)
d2vz1=V2_3(end)
d3vz1=V2_4(end)
vz2=V2_1(end)
dvz2=V2_2(end)
d2vz2=V2_3(end)
d3vz2=V2_4(end)
uz1
uz2
vz1
vz2
duz1
duz2
dvz1
dvz2
ux1=-duz1/(1i*k)
ux2=-duz2/(1i*k)
vx1=-dvz1/(1i*k)
vx2=-dvz2/(1i*k)
dux1=-d2uz1/(1i*k)
dux2=-d2uz2/(1i*k)
dvx1=-d2vz1/(1i*k)
dvx2=-d2vz2/(1i*k)
d2ux1=-d3uz1/(1i*k)
d2ux2=-d3uz2/(1i*k)
d2vx1=-d3vz1/(1i*k)
d2vx2=-d3vz2/(1i*k)
p1= (-re*vz1 + (-(a+g*i*k*z)*(re*(vx1/g)) + d2vx1 - k^2*vx1))/(i*k) % fluid pressure (from 1st set)
%dp1= (-(a+l*i*k*zb)*(re/l)+d2vz1-k^2*vz1)
pg1= (-(re/g)*a^2*ux1 + (1+er*a)*(d2ux1 - k^2*ux1))/(i*k) % Gel pressure (from 1st set)
%dpg1=(-(re/l)*l^2*uz1+(1+er*a)*(d2uz1-k^2*uz1))
p2= (-re*vz2 + (-(a+g*i*k*z)*(re*(vx2/g)) + d2vx2 - k^2*vx2))/(i*k) % fluid pressure (from 2nd set)
%dp2= (-(a+l*i*k*zb)*(re/l)+d2vz2-k^2*vz2)
pg2=(-(re/g)*a^2*ux1 + (1+er*a)*(d2ux2 - k^2*ux2))/(i*k) % Gel pressure (from 2nd set)
%dpg2=(-(re/l)*l^2*uz2+(1+er*a)*(d2uz2-k^2*uz2))
sigma_xz1= (1+er*a)*(dux1 + i*k*uz1) % SHEAR STRESS FOR THE GEL 1ST SET
sigma_xz2= (1+er*a)*(dux2 + i*k*uz2) % SHEAR STRESS FOR THE GEL 2ND SET
sigma_zz1= -pg1 + 2*(1+er*a)*duz1 % NORMAL STRESS FOR THE GEL 1ST SET
sigma_zz2= -pg2 + 2*(1+er*a)*duz2 % NORMAL STRESS FOR THE GEL 2ND SET
tou_xz1= dvx1 + i*k*vz1; % SHEAR STRESS FOR THE FLUID 1ST SET
tou_xz2= dvx2 + i*k*vz2; % SHEAR STRESS FOR THE FLUID 2ND SET
tou_zz1= -p1 + 2*dvz1; % NORMAL STRESS FOR THE FLUID 1ST SET
tou_zz2= -p2 + 2*dvz2; % NORMAL STRESS FOR THE FLUID 2ND SET
%% Declaring the elements of 4X4 MATRIX
a11= a*uz1
a12= a*uz2
a13= vz1
a14= vz2
a21= a*ux2 -g*uz2
a22= a*ux2 -g*uz2
a23= vx1
a24= vx2
a31= sigma_xz1
a32= sigma_xz2
a33= tou_xz1
a34= tou_xz2
a41= sigma_zz1
a42= sigma_zz2
a43= tou_zz1
a44= tou_zz2
A=[a11 a12 a13 a14; a21 a22 a23 a24; a31 a32 a33 a34; a41 a42 a43 a44]
F=det(A)% taking the determinant and setting it to the tolerance limit in Newton Raphson mentioned above
end
  4 个评论
显示 2更早的评论
SATYAJIT BHATTACHARJEE
SATYAJIT BHATTACHARJEE 2020-9-30
The code is as below
function newtonraphson
a0= -0.5 %input('enter initial guess for complex velocity of perturbations(c)')
itermax= 30 %input('enter max. iteration')
h=0.1 %input('enter value for stepsize(h)')
tolerance= 10^-10 %input('enter the tolerance limit')
glist=linspace(1,10,itermax); % Enter the values of g
for m=1:length(glist)
g=glist(m)
for j=1:itermax
F=f(a0,g)
Fprime=(f(a0+h,g)-f(a0,g))/h
a(j)=a0-(F/Fprime)
Fsolution=f(a(j),g)
if abs(Fsolution)<=tolerance
a(j)
break
else
a0=a(j)
end
end
alist(m)=a(end)
end
figure(1); plot(glist, real(alist), 'b')
xlabel('gamma')
ylabel('a-real')
title('plot of a-real vs gamma')
figure(2); plot(glist, imag(alist), 'r')
xlabel('gamma')
ylabel('a-imaginay')
title('plot of a-imaginary vs gamma')
end
function F=f(a,g)
syms uz(z) vz(z) Y
duz=diff(uz);
d2uz=diff(uz,2);
d3uz=diff(uz,3);
d4uz=diff(uz,4);
dvz=diff(vz);
d2vz=diff(vz,2);
d3vz=diff(vz,3);
d4vz=diff(vz,4);
er=0;% relative viscosity eta_r
H=0.1;% Ratio of thickness of gel to fluid
k=1;% gemma (dimensionless velocity)
re=0;% reynolds number
eqn1= (-k^2*d2uz + k^4*uz + i*k^3*a*er*d2uz - i*k^5*a*er*uz + g^2*k^4*uz - 2*i*k^3*g*duz - (re/g)*k^4*a^2*uz + d4uz - k^2*d2uz - i*k*a*er*d4uz + i*k^3*a*er*d2uz - g^2*k^2*d2uz + 2*i*k*g*d3uz + (re/g)*k^2*a^2*d2uz) ==0;
[UF] = odeToVectorField(eqn1);
coupled1 = matlabFunction(UF, 'Vars',{z,Y});
output1=coupled1;
uz10=0;
duz10=0;
d2uz10=1;
d3uz10=0;
zspan1=linspace(-H,0,10); % range of z
z10 = [uz10,duz10,d2uz10,d3uz10]; % array of initial values
[r1,v1] = ode45(@(z,Y)coupled1(z,Y),zspan1,z10)
%% Declaring the 2nd set independent solutions for displacement field at z=-10
uz20=0;
duz20=0;
d2uz20=0;
d3uz20=1;
zspan2=linspace(-H,0,10); % range of z
z20 = [uz20,duz20,d2uz20,d3uz20]; % array of initial values
[r2,v2] = ode45(@(z,Y)coupled1(z,Y),zspan2,z20)
eqn2= d4vz -2*k^2*d2vz + k^4*vz -i*re*k*z*d2vz - i*k*re*dvz + i*k*(re/g)*a*d2vz - re*d2vz + i*k^3*re*z*vz - (re/g)*i*k^3*a*vz;
[VF] = odeToVectorField(eqn2);
coupled2 = matlabFunction(VF, 'Vars',{z,Y});
output2=coupled2;
vz10=0;
dvz10=0;
d2vz10=1;
d3vz10=0;
Zspan1=linspace(1,0,10); % range of z
Z10 = [vz10,dvz10,d2vz10,d3vz10]; % array of initial values
[R1,V1] = ode45(@(z,Y)coupled2(z,Y),Zspan1,Z10)
%% Declaring the 2nd set independent solutions for velocity field at z=1
vz20=0;
dvz20=0;
d2vz20=0;
d3vz20=1;
Zspan2=linspace(1,0,10); % range of z
Z20 = [vz20,dvz20,d2vz20,d3vz20]; % array of initial values
[R2,V2] = ode45(@(z,Y)coupled2(z,Y),Zspan2,Z20)
v1_1=v1(:,1)
v1_2=v1(:,2)
v1_3=v1(:,3)
v1_4=v1(:,4)
v2_1=v2(:,1)
v2_2=v2(:,2)
v2_3=v2(:,3)
v2_4=v2(:,4)
V1_1=V1(:,1)
V1_2=V1(:,2)
V1_3=V1(:,3)
V1_4=V1(:,4)
V2_1=V2(:,1)
V2_2=V2(:,2)
V2_3=V2(:,3)
V2_4=V2(:,4)
uz1=v1_1(end)
duz1=v1_2(end)
d2uz1=v1_3(end)
d3uz1=v1_4(end)
uz2=v2_1(end)
duz2=v2_2(end)
d2uz2=v2_3(end)
d3uz2=v2_4(end)
vz1=V2_1(end)
dvz1=V2_2(end)
d2vz1=V2_3(end)
d3vz1=V2_4(end)
vz2=V2_1(end)
dvz2=V2_2(end)
d2vz2=V2_3(end)
d3vz2=V2_4(end)
% obtaining Values of variables at z=0
uz1
uz2
vz1
vz2
duz1
duz2
dvz1
dvz2
ux1=-duz1/(1i*k)
ux2=-duz2/(1i*k)
vx1=-dvz1/(1i*k)
vx2=-dvz2/(1i*k)
dux1=-d2uz1/(1i*k)
dux2=-d2uz2/(1i*k)
dvx1=-d2vz1/(1i*k)
dvx2=-d2vz2/(1i*k)
d2ux1=-d3uz1/(1i*k)
d2ux2=-d3uz2/(1i*k)
d2vx1=-d3vz1/(1i*k)
d2vx2=-d3vz2/(1i*k)
p1= (-(re*(vx1/g))*(g*i*k*z - a*k*i + g) + d2vx1 - k^2*vx1)/(i*k) % fluid pressure (from 1st set)
%dp1= (-(a+l*i*k*zb)*(re/l)+d2vz1-k^2*vz1)
pg1= ((1-i*k*a*er)*(d2ux1 - k^2*ux1) - g^2*k^2*ux1 + 2*i*k*g*dux1 + (re/g)*k^2*a^2*ux1)/(i*k) % Gel pressure (from 1st set)
%dpg1=(-(re/l)*l^2*uz1+(1+er*a)*(d2uz1-k^2*uz1))
p2= (-(re*(vx2/g))*(g*i*k*z - a*k*i + g) + d2vx2 - k^2*vx2)/(i*k) % fluid pressure (from 2nd set)
%dp2= (-(a+l*i*k*zb)*(re/l)+d2vz2-k^2*vz2)
pg2= ((1-i*k*a*er)*(d2ux2 - k^2*ux2) - g^2*k^2*ux2 + 2*i*k*g*dux2 + (re/g)*k^2*a^2*ux2)/(i*k) % Gel pressure (from 2nd set)
%dpg2=(-(re/l)*l^2*uz2+(1+er*a)*(d2uz2-k^2*uz2))
sigma_xz1= (1-i*k*a*er)*(dux1 + i*k*uz1) % SHEAR STRESS FOR THE GEL 1ST SET
sigma_xz2= (1-i*k*a*er)*(dux2 + i*k*uz2) % SHEAR STRESS FOR THE GEL 2ND SET
sigma_zz1= pg1 - 2*(1-i*k*er*a)*duz1 - 2*i*k*g*uz1 % NORMAL STRESS FOR THE GEL 1ST SET
sigma_zz2= pg2 - 2*(1-i*k*er*a)*duz2 - 2*i*k*g*uz2 % NORMAL STRESS FOR THE GEL 2ND SET
tou_xz1= dvx1 + i*k*vz1; % SHEAR STRESS FOR THE FLUID 1ST SET
tou_xz2= dvx2 + i*k*vz2; % SHEAR STRESS FOR THE FLUID 2ND SET
tou_zz1= p1 - 2*dvz1; % NORMAL STRESS FOR THE FLUID 1ST SET
tou_zz2= p2 - 2*dvz2; % NORMAL STRESS FOR THE FLUID 2ND SET
%% Declaring the elements of 4X4 MATRIX
a11= -i*k*a*uz1
a12= -i*k*a*uz2
a13= vz1
a14= vz2
a21= -i*k*a*ux1 -g*uz1
a22= -i*k*a*ux2 -g*uz2
a23= vx1
a24= vx2
a31= sigma_xz1
a32= sigma_xz2
a33= tou_xz1
a34= tou_xz2
a41= sigma_zz1
a42= sigma_zz2
a43= tou_zz1
a44= tou_zz2
A=[a11 a12 a13 a14; a21 a22 a23 a24; a31 a32 a33 a34; a41 a42 a43 a44]
F=det(A)% taking the determinant and setting it to the tolerance limit in Newton Raphson mentioned above
end
VBBV
VBBV 2020-9-30
编辑：VBBV 2020-9-30
Use the break in the else part of if else end
%if true
% code
% end
if condition
...
else
...
break;
end
Also the function called F =f(a,g) within function newtonraphson, Check the size of zspan it is divided in 10 intervals ... set it to itermax

请先登录，再进行评论。

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

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Stress and Strain 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by