Why i don not get the graph when i add A1(x) in the function f1=@(x,y) A1(x)+(1.3​3*pi*eta*y​^3)/(1+1.3​3*pi*eta*y​^3);

1 view (last 30 days)
%clc
%clear all
eta=0.25;
L=2.5;
a=0.1;
h=0.1;
b=6;
x=a:h:b;
%%%%%%%%%%%%%%%%%%
A =@(x) (0<x).*1+(0<=x & x<=L).*(1+0.5*(1-cos(2*pi*x/L)).^1/2) + (L<x).*1;
A1=@(x) (gradient(A(x))./gradient(x));
f1=@(x,y) A1(x)+(1.33*pi*eta*y^3)/(1+1.33*pi*eta*y^3);
y(1)=01;
for i=1:length(x)-1;
k1=h*f1(x(i),y(i));
k2=h*f1(x(i)+0.5*h,y(i)+.5*k1);
k3=h*f1(x(i)+.5*h,y(i)+.5*k2);
k4=h*f1(x(i)+h,y(i)+k3);
y(i+1)=y(i)+(1/6)*(k1+2*k2+2*k3+k4);
end
plot(x,y)
  4 Comments
Farooq Hussain
Farooq Hussain on 23 Apr 2020
Thank you again very much for this latest suggestion...But when i apply the same procedure i get the same plot whether i include or remove A(x) in f3 of the code....Can you please explain
clc
clear all;
%%%%
L=1.0;
Re=1000;
k=1.4;
we=137;
d=0.8; %%%% d=sigma
alphas=0;
%%% Step size%%%
a1=0;
h=0.05;
n=5;
x=[a1:h:n];
%%% Functions %%%
A = @(x)(0<x).*1+(0<=x & x<=L).*(1+0.5*(1-cos(2*pi*x/L))).^1/2 + (L<x).*1;
f1=@(x,u,v,cp,H,R)v;
f2=@(x,u,v,cp,H,R)H;
f3=@(x,u,v,cp,H,R)(-1+alphas+alphas*R.^3)+A(x).*u;
f4=@(x,u,v,cp,H,R)(2*u*v*(1-alphas))/(-1+alphas+alphas*R^3);
f5=@(x,u,v,cp,H,R)(2*(R^(-3*k)-(R^-1))/we+0.5*d*(-1+R^(-3*k))-0.5*cp-u*v*R*H-1.5*(u^2)*(H^2)-(4*u*R*H)/(Re))./(R*u^2);
%%% Initial Conditions %%%
cp(1)=0;
u(1)=01;
v(1)=0.0154;
R(1)=01;
H(1)=0.021;
for i=1:length(x)-1;
m1=h*f1(x(i), u(i), v(i), cp(i), H(i), R(i));
n1=h*f2(x(i), u(i), v(i), cp(i), H(i), R(i));
p1=h*f3(x(i), u(i), v(i), cp(i), H(i), R(i));
q1=h*f4(x(i), u(i), v(i), cp(i), H(i), R(i));
r1=h*f5(x(i), u(i), v(i), cp(i), H(i), R(i));
m2=h*f1(x(i)+0.5*h, u(i)+0.5*m1, v(i)+0.5*n1, cp(i)+0.5*p1, H(i)+0.5*q1, R(i)+0.5*r1);
n2=h*f2(x(i)+0.5*h, u(i)+0.5*m1, v(i)+0.5*n1, cp(i)+0.5*p1, H(i)+0.5*q1, R(i)+0.5*r1);
p2=h*f3(x(i)+0.5*h, u(i)+0.5*m1, v(i)+0.5*n1, cp(i)+0.5*p1, H(i)+0.5*q1, R(i)+0.5*r1);
q2=h*f4(x(i)+0.5*h, u(i)+0.5*m1, v(i)+0.5*n1, cp(i)+0.5*p1, H(i)+0.5*q1, R(i)+0.5*r1);
r2=h*f5(x(i)+0.5*h, u(i)+0.5*m1, v(i)+0.5*n1, cp(i)+0.5*p1, H(i)+0.5*q1, R(i)+0.5*r1);
m3=h*f1(x(i)+0.5*h, u(i)+0.5*m2, v(i)+0.5*n2, cp(i)+0.5*p2, H(i)+0.5*q2, R(i)+0.5*r2);
n3=h*f2(x(i)+0.5*h, u(i)+0.5*m2, v(i)+0.5*n2, cp(i)+0.5*p2, H(i)+0.5*q2, R(i)+0.5*r2);
p3=h*f3(x(i)+0.5*h, u(i)+0.5*m2, v(i)+0.5*n2, cp(i)+0.5*p2, H(i)+0.5*q2, R(i)+0.5*r2);
q3=h*f4(x(i)+0.5*h, u(i)+0.5*m2, v(i)+0.5*n2, cp(i)+0.5*p2, H(i)+0.5*q2, R(i)+0.5*r2);
r3=h*f5(x(i)+0.5*h, u(i)+0.5*m2, v(i)+0.5*n2, cp(i)+0.5*p2, H(i)+0.5*q2, R(i)+0.5*r2);
m4=h*f1(x(i)+h, u(i)+m3, v(i)+n3, cp(i)+p3, H(i)+q3, R(i)+r3);
n4=h*f2(x(i)+h, u(i)+m3, v(i)+n3, cp(i)+p3, H(i)+q3, R(i)+q3);
p4=h*f3(x(i)+h, u(i)+m3, v(i)+n3, cp(i)+p3, H(i)+q3, R(i)+r3);
q4=h*f4(x(i)+h, u(i)+m3, v(i)+n3, cp(i)+p3, H(i)+q3, R(i)+r3);
r4=h*f5(x(i)+h, u(i)+m3, v(i)+n3, cp(i)+p3, H(i)+q3, R(i)+r3);
u(i+1)=u(i)+(1/6)*(m1+2*m2+2*m3+m4);
v(i+1)=v(i)+(1/6)*(n1+2*n2+2*n3+n4);
cp(i+1)=cp(i)+(1/6)*(p1+2*p2+2*p3+p4);
H(i+1)=H(i)+(1/6)*(q1+2*q2+2*q3+q4);
R(i+1)=R(i)+(1/6)*(r1+2*r2+2*r3+r4);
end
plot(x,u,'k')
Star Strider
Star Strider on 23 Apr 2020
The ‘f3’ function likely makes a very small contribution, regardless. (The code is so complicated that it is difficult to follow.)
Plotting ‘f3’ alone (after the end of the loop so that all the vectors are present) with:
f3=@(x,u,v,cp,H,R)(-1+alphas+alphas*R.^3)+A(x).*u;
or without:
f3=@(x,u,v,cp,H,R)(-1+alphas+alphas*R.^3);%+A(x).*u;
the ‘A(x)’ term definitely makes a difference:
f3v = f3(x+h, u+m3, v+n3, cp+p3, H+q3, R+r3);
figure
plot(x, h*f3v);
So ‘f3’ works correctly.
.

Sign in to comment.

Accepted Answer

Star Strider
Star Strider on 18 Apr 2020
Do not use the for loop here. Let vectorisation work for you.
Try this:
eta=0.25;
L=2.5;
a=0.1;
h=0.1;
b=6;
x=a:h:b;
%%%%%%%%%%%%%%%%%%
A =@(x) (0<x).*1+(0<=x & x<=L).*(1+0.5*(1-cos(2*pi*x/L)).^1/2) + (L<x).*1;
A1=@(x) (gradient(A(x))./gradient(x));
f1=@(x,y) A1(x)+(1.33*pi*eta*y.^3)./(1+1.33*pi*eta*y.^3);
y(1)=01;
% for i=1:length(x)-1;
k1=h*f1(x,y);
k2=h*f1(x+0.5*h,y+.5*k1);
k3=h*f1(x+.5*h,y+.5*k2);
k4=h*f1(x+h,y+k3);
y = +(1/6)*(k1+2*k2+2*k3+k4);
% end
plot(x,y)
.

More Answers (0)

Tags

Community Treasure Hunt

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

Start Hunting!