MATLAB Answers

symbolic and numerical code merged to solve

6 views (last 30 days)
function main
%%%%%%NUMERICAL CODE
A=0.5; pr=1; a=1; phi=0.1;a1=2;a2=1;xa=0;xb=6;
solinit=bvpinit(linspace(xa,xb,10),[0 1 0 a 1 0 0 0]);
sol=bvp5c(@ode,@bc,solinit);
x=linspace(xa,xb,100);S=deval(sol,x);
function res=bc(ya,yb)
res=[ya(1); ya(2)-1; ya(4); ya(5)-a; ya(7)-1; yb(2); yb(5); yb(7)];
end
function dydx=ode(x,y)
dydx=[y(2); y(3); 2*a1*y(2)*(y(2)+y(5))-a1*y(3)*(y(1)+y(4));
y(5); y(6); 2*a1*y(5)*(y(2)+y(5))-a1*y(6)*(y(1)+y(4));
y(8); A*pr*a2*y(7)*(y(2)+y(5))-pr*a2*y(8)*(y(1)+y(4))];
end
f0 = deval(sol,0);
p=f0(3);q=f0(6);r=f0(8);
figure(1)
plot(x,S([2],:)); %for f'
xlabel('\eta'); ylabel('f`');
hold on
% %%%%%%%%%%%%%%%%% SYMBOLIC CODE
syms t x a p q r a1 a2 A pr
f(1)=x+(p/2)*x^2;g(1)=a*x+(q/2)*x^2;h(1)=1+r*x;
for i=1:3
fa(i) = subs(f(i),x,t);dfa = diff(fa(i),t,1);d2fa = diff(dfa,t,1);
ga(i) = subs(g(i),x,t);dga = diff(ga(i),t,1);d2ga = diff(dga,t,1);
ha(i) = subs(h(i),x,t);dha = diff(ha(i),t,1);
If1=int((fa(i)+ga(i))*(-d2fa) +2*dfa*(dfa+dga),t,0,t);If2=int(If1,t,0,t);If3=int(If2,t,0,x);
Ig1=int((fa(i)+ga(i))*(-d2ga) + 2*dfa*(dfa+dga),t,0,t);Ig2=int(Ig1,t,0,t);Ig3=int(Ig2,t,0,x);
Ih1=int((fa(i)+ga(i))*(-dha) + A*ha(i)*(dfa+dga),t,0,t);Ih2=int(Ih1,t,0,x);
f(i+1) = a1*If3;g(i+1) = a1*Ig3;h(i+1) = pr*a2*Ih2;
disp(f(i+1))
end
f=f(1)+f(2)+f(3);g=g(1)+g(2)+g(3);h=h(1)+h(2)+h(3);
x=0.0:0.01:6;
F=[0 diff(f)];
figure(2)
fplot(x,F,'LineWidth',1.5)
hold on
end
%%%%%%%%%
I want to run the symbolic code taking the values of p, q, r from NUMERIC CODE and also to draw fig(2) but gives error. Also to check the expression f=f(1)+f(2)+f(3); is right or not.
Any help will be appreciated.

  0 Comments

Sign in to comment.

Accepted Answer

Walter Roberson
Walter Roberson on 22 Dec 2019
Guessing about what you want:
function main
%%%%%%NUMERICAL CODE
A=0.5; pr=1; a=1; a1=2;a2=1;xa=0;xb=6;
solinit=bvpinit(linspace(xa,xb,10),[0 1 0 a 1 0 0 0]);
sol=bvp5c(@ode,@bc,solinit);
xn=linspace(xa,xb,100);
x = xn;
S=deval(sol,x);
function res=bc(ya,yb)
res=[ya(1); ya(2)-1; ya(4); ya(5)-a; ya(7)-1; yb(2); yb(5); yb(7)];
end
function dydx=ode(~,y)
dydx=[y(2); y(3); 2*a1*y(2)*(y(2)+y(5))-a1*y(3)*(y(1)+y(4));
y(5); y(6); 2*a1*y(5)*(y(2)+y(5))-a1*y(6)*(y(1)+y(4));
y(8); A*pr*a2*y(7)*(y(2)+y(5))-pr*a2*y(8)*(y(1)+y(4))];
end
f0 = deval(sol,0);
p=f0(3);q=f0(6);r=f0(8);
figure(1)
plot(x,S([2],:)); %for f'
xlabel('\eta'); ylabel('f`');
hold on
% %%%%%%%%%%%%%%%%% SYMBOLIC CODE
t = sym('t');
x = sym('x');
%{
p = sym('p');
q = sym('q');
r = sym('r');
a1 = sym('a1');
a2 = sym('a2');
A = sym('A');
pr = sym('pr');
%}
f = zeros(1,3,'sym');
g = zeros(1,3,'sym');
h = zeros(1,3,'sym');
fa = zeros(1,3,'sym');
ga = zeros(1,3,'sym');
ha = zeros(1,3,'sym');
f(1)=x+(p/2)*x^2;g(1)=a*x+(q/2)*x^2;h(1)=1+r*x;
for i=1:3
fa(i) = subs(f(i),x,t);dfa = diff(fa(i),t,1);d2fa = diff(dfa,t,1);
ga(i) = subs(g(i),x,t);dga = diff(ga(i),t,1);d2ga = diff(dga,t,1);
ha(i) = subs(h(i),x,t);dha = diff(ha(i),t,1);
If1=int((fa(i)+ga(i))*(-d2fa) +2*dfa*(dfa+dga),t,0,t);If2=int(If1,t,0,t);If3=int(If2,t,0,x);
Ig1=int((fa(i)+ga(i))*(-d2ga) + 2*dfa*(dfa+dga),t,0,t);Ig2=int(Ig1,t,0,t);Ig3=int(Ig2,t,0,x);
Ih1=int((fa(i)+ga(i))*(-dha) + A*ha(i)*(dfa+dga),t,0,t);Ih2=int(Ih1,t,0,x);
f(i+1) = a1*If3;g(i+1) = a1*Ig3;h(i+1) = pr*a2*Ih2;
disp(f(i+1))
end
f=f(1)+f(2)+f(3);g=g(1)+g(2)+g(3);h=h(1)+h(2)+h(3);
F=[0 diff(double(subs(f,x,xn)))];
figure(2)
plot(xn,F,'LineWidth',1.5)
hold on
end
%%%%%%%%%

  8 Comments

Show 5 older comments
MINATI
MINATI on 26 Dec 2019
But if I put '1' in that place, same error happens
MINATI
MINATI on 26 Dec 2019
And for diff(f), any changes.
The code for 'h' also didnot obey B.C(ya(7)-1;).

Sign in to comment.

More Answers (1)

ikram
ikram on 18 Dec 2019
syms t x a p q r a1 a2 A pr
you define variable use below code
t=sym('t')
syms x a p q r a1 a2 A pr
but your code still can't run correctly

  10 Comments

Show 7 older comments
MINATI
MINATI on 31 Dec 2019
No it is not matching, I had checked earlier
Walter Roberson
Walter Roberson on 31 Dec 2019
It works when I try it.
xlim([0 1])
497235.jpg
It is pretty clear that F is starting from 1, as required.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!