symbolic and numerical code merged to solve

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.

Answer 1

[Walter Roberson]

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
%%%%%%%%%

[MINATI]

@Dear Walter

Both the figs should start from '1' initially (As ya(2)-1;) but fig. 2 did not obey.

[Walter Roberson]

No, the figures should not start at 1 initially. You use linspace to create x between xa=0 and xb=6 for your first figure does plot(x,S([2],:)); so it has to start from 0. For your second figure, you use x=0.0:0.01:6 which is not necessarily exactly the same as the linspace() but again starts from 0, not from 1.

If you are referring to the y value, then look again at your code:

F=[0 diff(f)];
figure(2)
fplot(x,F,'LineWidth',1.5)

You hard-coded a 0 at the beginning of F, so the plot cannot start from 1.

[MINATI]

@Dear Walter

diff(f)=1 at x=0 which obeys B.C ( ya(2)-1) but fig.2 shows another behaviour. I am not able to catch the mistake.

And also to draw (h), any modifications needed?

[Walter Roberson]

If f'(0) is 1 then why do you hardcode 0 at the beginning of F?

[Walter Roberson]

 H = double(subs(h, x, xn)) ;
 hold on
 plot(xn, H)

[MINATI]

But if I put '1' in that place, same error happens

[MINATI]

And for diff(f), any changes.

The code for 'h' also didnot obey B.C(ya(7)-1;).

[MINATI]

ok

I will try the rest.

Thanks

Answer 2

[ikram]

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

[MINATI]

Both the codes individualiy runs correctly but I want to retrieve the values of p,q,r from numerical code to symbolic code.

Any idea

[Walter Roberson]

When you have

syms a

That is the same as

a = sym('a');

which replaces the current a with a reference to a variable a that lives inside the symbolic engine. This overwrites the a you had with a=1; in the numeric code, leaving a as an undefined variable. That variable shows up as an unresolved variable in the expression you are trying to fplot()

[MINATI]

since f '(0)=1; so 'F' should start from '1' in Y-axis but here it is not, why?

[Walter Roberson]

You hard-code

    F=[0 diff(double(subs(f,x,xn)))];

so F will start with 0.

[MINATI]

F=[0 diff(double(subs(f,x,xn)))];

Can't that zero(0) be changed to '1' OR any other number, say, 'a'

[Walter Roberson]

Yes

[MINATI PATRA]

How Because in purely symbolic code by replacing 1, the curve starts from 1 but here it is not.

[Walter Roberson]

   F=[a diff(double(subs(f,x,xn)))];

[MINATI]

No it is not matching, I had checked earlier

[Walter Roberson]

It works when I try it.

xlim([0 1])

It is pretty clear that F is starting from 1, as required.