Items tagged with code code Tagged Items Feed

Dear all,

I would like to solve the Fredholm Integral equation, using numerical method.
This is my code.

there is a problem with subs, does not working here.
# Then, we obtain from (9) the coeficient A[n] and B[n].

Then I woulk like to recompute (2), and then compute (1).
# Puting x=m*h, in (1), how can we generate a linear Matrix from (1).

 

Fred.mw

Thanks

Dear all;

Please in this code Maple, I plot using pointplot3d, some points. I would like to add a name of each point.

The name of each point  P:=[x(i),y(j),t(k)] is u[i,j]^k.  Please how can I add the name at each point. here is the maple code. 

Thanks.

pointplot3d.mw

Dear all;

 

I need a help in this question.

u(x,y,t) my solution of PDE. x,y space, and t=time.

In the case without t. .i.e. u(x,y). Here is a visualization of the lattice u(x_i,y_j). i=1..3, anf j=1..3.  Please try this example, it's working. 

my question, if i would like to add a third variable t, i.e.  and get u(x_i,y_j,t_k)  on each points. I would like to show the lattice in (x,y) plan for each t_k. I Think I will get many parallel lattice. 

Please can some one, modify this code to get the parallel lattice. Thanks.


with(LinearAlgebra):
with(plots):
with(PDEtools):
interface(rtablesize=20):

Here,


L := 'L':
N := 'N':
g := 'g':
Z := i -> -L+2*L/(N+1)*i;
x[0] = Z(0),x[N+1] = Z(N+1),y[0] = Z(0),y[N+1] = Z(N+1);

 

N := 4;
L := 1;
r := L/(N+1)/4;
ngon := (n,x,y,r,phi) -> [seq([x+r*cos(2*Pi*i/n+phi), y+r*sin(2*Pi*i/n+phi)], i = 1 .. n)]:
p[1] := display([seq(polygonplot(ngon(4,Z(0),Z(j),r,Pi/2),color=magenta),j=0..N+1),

seq(polygonplot(ngon(4,Z(N+1),Z(j),r,Pi/2),color=magenta),j=0..N+1),

seq(polygonplot(ngon(4,Z(i),Z(0),r,0),color=magenta),i=1..N),

seq(polygonplot(ngon(4,Z(i),Z(N+1),r,0),color=magenta),i=1..N),

seq(seq(polygonplot(ngon(20,Z(i),Z(j),r,0),color=white),i=1..N),j=1..N),

textplot([seq(seq([Z(i+0.1),Z(j),typeset(u[i,j])],i=0..N+1),j=0..N+1)],align={above,right})],view=[Z(-1)..Z(N+2),Z(-1)..Z(N+2)],tickmarks=[[seq(Z(i)=typeset(x[i]=evalf[2](Z(i))),i=0..N+1)],[seq(Z(i)=typeset(y[i]=evalf[2](Z(i))),i=0..N+1)]],

axes=boxed,scaling=constrained,labels=[``,``]):
p[1];

 

 

 

 

 

 

Hi all

I dont know why some Z1 appears on the screen and the code does not converge.
please help me
thanks alooooot


restart;
n:=3;
nn:=3;
m:=1;
BB:=1;
BINF:=5:
pr:=7;
digits:=10;
>
eq1:=diff(f(tau),tau$3)+((3/5)*f(tau)*diff(f(tau),tau$2))-(1/5)*(diff(f(tau),tau$1))^2+((2/5)*tau*diff(h(tau),tau$1))-((2/5)*h(tau))-BB*diff(f(tau),tau$1)=0;
eq11:=(1/pr)*diff(h(tau),tau$3)+(3/5)*f(tau)*diff(h(tau),tau$2)=0;

h(tau):=sum(p^i*h[i](tau),i=0..nn);
f(tau):=sum(p^i*f[i](tau),i=0..n);

 

H1:= p*(diff(f(tau),tau$3)+((3/5)*f(tau)*diff(f(tau),tau$2))-(1/5)*(diff(f(tau),tau$1))^2+((2/5)*tau*diff(h(tau),tau$1))-((2/5)*h(tau))-BB*diff(f(tau),tau$1))+(1-p)*(diff(f(tau),tau$3)):
H11:= p*((1/pr)*diff(h(tau),tau$3)+(3/5)*f(tau)*diff(h(tau),tau$2))+(1-p)*(diff(h(tau),tau$3)):
>
eq2:=simplify(H1):
eq22:=simplify(H11):
eq3:=collect(expand(eq2),p):
eq33:=collect(expand(eq22),p):
eq4:=
convert(series(collect(expand(eq2), p), p, n+1), 'polynom');
eq44:=
convert(series(collect(expand(eq22), p), p, n+1), 'polynom');
for i to n do
s[i] := coeff(eq4, p^i) ;
print (i);
end do:
for i to nn do
ss[i] := coeff(eq44, p^i) ;
print (i);
end do:
s[0]:=diff(f[0](tau), tau$3);
ss[0]:=diff(h[0](tau), tau$3);
ics[0]:=f[0](0)=0, D(f[0])(0)=0, D(f[0])(BINF)=0;
icss[0]:=h[0](BINF)=0, D(h[0])(0)=1, D(h[0])(BINF)=0;

dsolve({s[0], ics[0]});
f[0](tau):= rhs(%);
#dsolve({ss[0], icss[0]});
h[0](tau):= -exp(-tau); #;rhs(%);

>
>
for i to n do
f[ii-1](tau):=convert(series(f[ii-1](tau), tau, nn+1), 'polynom');
h[ii-1](tau):=convert(series(h[ii-1](tau), tau, nn+1), 'polynom');
s[i]:=simplify((s[i]));
ics[i]:=f[i](0)=0, D(f[i])(0)=0, D(f[i])(BINF)=0;
dsolve({s[i], ics[i]});
f[i](tau):=rhs(%);
ss[i]:=(ss[i]);
icss[i]:=h[i](BINF)=0, D(h[i])(0)=0, D(h[i])(BINF)=0;
dsolve({ss[i], icss[i]});
h[i](tau):=rhs(%);

end do;

f(tau):=sum((f[j])(tau),j=0..n);
with(numapprox):

 

plot(diff(f(tau),tau),tau=0..5,color=blue,style=point,symbol=circle,symbolsize=7,labels=["tau","velocity"]);
plot(pade(diff(f(tau),tau), tau, [7, 7]),tau=0..5,color=blue,style=point,symbol=circle,symbolsize=7,labels=["tau","velocity"]);

 

 

Hi.

please, I need some code about stepsize adaptative Runge Kutta method ( any order ) or other methods.

Thanks

Dear all;

Please I need your help to find the error in my code.

I want to solve an ode, with condition on step size.

ode := diff(y(x), x) = 2*x+y(x);
f:=(x,y)->2*x-y;

analyticsol := rhs(dsolve({ode, y(0) = 1}));
RKadaptivestepsize := proc (f, a, b, epsilon, N)
local x, y, n, h,k,z,R,p;
p:=2;
h := evalf(b-a)/N; ## we begin with this setpsize
x[0] := a; y[0] := 1; ## Initialisation
for n from 0 to N-1 do  ##loop
x[n+1] := a+(n+1)*h;  ## noeuds
k[1] := f(x[n], y[n]);
k[2] := f(x[n]+h, y[n]+h*k[1]);
k[3] := f(x[n]+h/2, y[n]+h/4*(k[1]+k[2]));
z[n+1] := z[n]+(h/2)*(k[1]+k[2]);## 2-stage runge Kutta.
y[n+1] := y[n]+(h/6)*(k[1]+k[2]+4*k[3]);
R:=abs(y[n+1]-z[n+1]); ## local erreur
hstar:=sqrt(epsilon/R)
if R=<=epsilon    then
   x[n] := x[n+1]+h;
   y[n]:=y[n+1];
   n:=n+1;

else

h:=hstar;
end if
 end do;
[seq([x[n], y[n]], n = 0 .. N)];
[seq([x[n], z[n]], n = 0 .. N)];
end proc:

epsilon:=1e-8;
RKadaptivestepsize((x,y)->2*x-y,0,1, epsilon,20)

Dear all

 

Please, I need you help, to make a simple procedure with two output.  I try this simple code :

test:= proc(a,b)

local variable .....;

for i from 1 to 10 do

expression

end do

retur  fff;

return ggg;

end proc

Then I do  : test(a,b) but I get only ggg and not fff

 

Any help please

 

Dear all,

Thank you for your Help.

h: stepsize;

x in [0,x0];

I give all the step of my code, but I think there is a mistake. I wait for your Help.

I would like to compute the error between  Method Huen with step size h and step size 2h using the definition of epsilon given below:

 ## The error written epsilon(x0,h)= sqrt(1/(N+1) * sum_i=0^N  (y_i^{2h}-y_(2i)^h)^2 ), where y_i^(2h) is the approximation of y(i*2*h).

 ## We want : loglog epsilon versus h.

  epsilon:=(x0,h)->sqrt( 1/(N+1)*add( (Heun1(f,x0,i)-Heun2(f,x0,i))^2,i=0..N ) );

  f:=(x,y)=1/(1+cos(y)); 

  ode:=diff(y(x),x)=f(x,y);

ic:=y(0)=1;  h:=x0/(2*N);

## Method Heun with step size 2h

> Heun1 := proc (f, x0,)

local x, y, i, h, k;

y := Array(0 .. N);

x := Array(0 .. N);

h := evalf((1/2)*x0/N);

x[0] := 0;

y[0] := 1;

for i from 0 to N do

x[i+1] := (2*i+2)*h;

k[1] := f(x[i], y[i]);

k[2] := f(x[i]+h, y[i]+h*k[1]);

y[i+1] := y[i]+h*((1/2)*k[1]+(1/2)*k[2]);

end do;

[seq([x[i], y[i]], i = 0 .. N)];

end proc;

### Now Heun with step size h  ( the same h)

> Heun2 := proc (f, x0,)

local x, y, i, h, k;

y := Array(0 .. N);

x := Array(0 .. N);

h := evalf((1/2)*x0/N);

x[0] := 0;

y[0] := 1;

for i from 0 to N do

x[i+1] := (i+1)*h;

k[1] := f(x[i], y[i]);

k[2] := f(x[i]+h, y[i]+h*k[1]);

y[i+1] := y[i]+h*((1/2)*k[1]+(1/2)*k[2]);

end do;

[seq([x[2*i], y[2*i]], i = 0 .. N)];

end proc;

 

 

Thanks you for your help.


                                

                        

 

i_before_e$.mw

My Homies

I'm trying to find words which contradict the rule "i before e except after c". according to BBC program "QI" there are 923 exceptions, but below i get less than that. i guess its only as good as the dictionary, the 'inbuilt' has only 25k words but 'ospd3' (with 80k words) seems to stop before letter Q (with my sledgehammer code) when there are more words out there like "weird", "veil" and "reciept".

with(StringTools); with(PatternDictionary); bid := Create(ospd3);
words :=([seq])(LowerCase(Get(bid, i)), i = 1 .. Size(bid)-1);
tot := 0;
for n to nops(words) do
C := searchtext(cie, words[n]); E := searchtext(ei, words[n]);
if`and`(C = 0, E = 0) then next else tot := tot+1;
print([tot, words[n]]) end if
end do

Hi,

I wrote the following code which is properly run

 


restart:

# parametrs

MUR:=(1-phi)^2.5:
RhoUR:=(1-phi+phi*rho[p]/rho[f]):
RhoCPR:=(1-phi+phi*rhocp[p]/rhocp[f]):
BetaUR:=(phi*rho[p]*beta[p]+(1-phi)*rho[f]*beta[f])/(RhoUR*rho[f])/beta[f]:

dqu3:=diff(h(x),x$1)-RhoUR*BetaUR*T(x);
dqu2:=5*diff(T(x),x$2)+k[f]/k[nf]*Pr*RhoCPR*f(x)*diff(T(x),x$1);
dqu1:=5/(MUR)*diff(f(x),x$3)
+ 2*(diff(h(x),x$1)*x-h(x))
+RhoUR*(3*f(x)*diff(f(x),x$2)-diff(f(x),x$1)^2);
rho[f]:=998.2: cp[f]:=4182: k[f]:=0.597:   beta[f]:= 2.066/10000:
rho[p]:=3380: cp[p]:=773: k[p]:=36:   beta[p]:= 8.4/1000000:

k[nf]:=((k[p]+2*k[f])-2*phi*(k[f]-k[p]))/((k[p]+2*k[f])+phi*(k[f]-k[p])):
rhocp[nf]:=rho[p]*cp[p]*phi+rho[f]*cp[f]*(1-phi):
rhocp[p]:=rho[p]*cp[p]:
rhocp[f]:=rho[f]*cp[f]:

phi:=0.00:
binfinitive:=6: Pr:=7: lambda:=0:


with(plots):
pppe:=dsolve( {dqu1=0,dqu2=0,dqu3=0,T(0)=1,T(binfinitive)=0,f(0)=0,D(f)(0)=lambda,D(f)(binfinitive)=0,h(binfinitive)=0}, numeric );
-pppe(0);
print(odeplot(pppe,[x,diff(f(x),x)],0..binfinitive,color=black,numpoints=400));
print(odeplot(pppe,[[x,diff(f(x),x)]],0..binfinitive,color=black,numpoints=400));
print(odeplot(pppe,[[x,T(x)]],0..binfinitive,color=black,numpoints=400));


However, in some range of parameters, I must increase the value of binfinitive (for example binfinitive=50). however, my code is doesnt converge for higher values of 10 (at most). Can anyone change this algorithm in a way that it insensitive to the value of binfinitive?

Many thanks for your attention in advance

 

Amir

Hi,

I have a problem with dsolve in the following code

restart;
>
n:=20;
m:=1;
cc:=-200;
zzeta:=0.1;
sefr1:=0.3;
sefr:=0.2;
MM:=0;
lambda:=0.1;
Br:=1;
nn:=3;
>
>
#u(tau):=tau;
u(tau):=421.7129935*tau-2217.587728*tau^2+8897.376593*tau^3-27612.59182*tau^4+64248.00336*tau^5-1.083977605*10^5*tau^6-10.57029600-1.080951714*10^6*tau^13+7.999517316*10^5*tau^14-4.788741005*10^5*tau^15+2.309563748*10^5*tau^16+26511.11102*tau^18-5959.001794*tau^19+1.148523882*10^5*tau^7-95.23809524*tau^21+4.545454545*tau^22-9435.563781*tau^8-2.587683745*10^5*tau^9+6.473880128*10^5*tau^10+948.0272727*tau^20-88660.41892*tau^17-1.008692404*10^6*tau^11+1.175504242*10^6*tau^12;
>
>
B := 1+(2*(1-zzeta))*Br*(int(tau*(diff(u(tau), tau))^2, tau = (1/2)*zzeta*(1-zzeta) .. 1/2*(1-zzeta)));
eq4 := 4*B*u(tau)-(1+zzeta)*(diff(tau*(diff(theta(tau), tau)), tau))/tau-(1+zzeta)*Br*(diff(u(tau), tau))^2;


theta(tau):=sum(p^ii*theta[ii](tau),ii=0..nn);
HH:= p*((4*(1+(2*(1-zzeta))*Br*(int(tau*(diff(u(tau), tau))^2, tau = (1/2)*zzeta*(1-zzeta) .. 1/2-(1/2)*zzeta))))*u(tau)-(1+zzeta)*(diff(theta(tau), tau)+tau*(diff(theta(tau), tau, tau)))/tau-(1+zzeta)*Br*(diff(u(tau), tau))^2)+(1-p)*(diff(theta(tau),tau$2)):
eq5:=simplify(HH):
eq6:=collect(expand(eq5),p);

eq7:=
convert(series(collect(expand(eq5), p), p, nn+1), 'polynom');


for ii to nn do
ss[ii] := (coeff(eq7, p^ii)) ;
print (ii);
end do;

ss[0]:=diff(theta[0](tau), tau, tau);

icss[0]:=theta[0](zzeta/(2*(1-zzeta)))=0, D(theta[0])(1/(2*(1-zzeta)))=1;

dsolve({ss[0], icss[0]});
theta[0](tau):= rhs(%);


for ii to nn do
ss[ii]:=evalf[5](ss[ii]);
icss[ii]:=theta[ii](zzeta/(2*(1-zzeta)))=0, D(theta[ii])(1/(2*(1-zzeta)))=0;
dsolve({ss[ii], icss[ii]});
theta[ii](tau):=rhs(%);
end do;

I would be most grateful if you help me to find this problem.

Thanks for your attention in advance

 

Dear all

I would like to convert Matlab code to Maple, is there anu idea, this is the code.

 

% Usage: [y t] = abm4(f,a,b,ya,n) or y = abm4(f,a,b,ya,n)
% Adams-Bashforth-Moulton 4-th order predictor-corrector method for initial value problems
% It uses
% Adams-Bashforth 4-step method as a precdictor,
% Adams-Moulton 3-step method as a corrector, and
% Runge-Kutta method of order 4 as a starter
%
% Input:
% f - Matlab inline function f(t,y)
% a,b - interval
% ya - initial condition
% n - number of subintervals (panels)
%
% Output:
% y - computed solution
% t - time steps
%
% Examples:
% [y t]=abm4(@myfunc,0,1,1,10);          here 'myfunc' is a user-defined function in M-file
% y=abm4(inline('sin(y*t)','t','y'),0,1,1,10);
% f=inline('sin(y(1))-cos(y(2))','t','y');
% y=abm4(f,0,1,1,10);

function [y t] = abm4(f,a,b,ya,n)
h = (b - a) / n;
h24 = h / 24;

y(1,:) = ya;
t(1) = a;

m = min(3,n);

for i = 1 : m % start-up phase, using Runge-Kutta of order 4
    t(i+1) = t(i) + h;
    s(i,:) = f(t(i), y(i,:));
    s2 = f(t(i) + h / 2, y(i,:) + s(i,:) * h /2);
    s3 = f(t(i) + h / 2, y(i,:) + s2 * h /2);
    s4 = f(t(i+1), y(i,:) + s3 * h);
    y(i+1,:) = y(i,:) + (s(i,:) + s2+s2 + s3+s3 + s4) * h / 6;
end;

for i = m + 1 : n % main phase
    s(i,:) = f(t(i), y(i,:));
    y(i+1,:) = y(i,:) + (55 * s(i,:) - 59 * s(i-1,:) + 37 * s(i-2,:) - 9 * s(i-3,:)) * h24; % predictor
    t(i+1) = t(i) + h;
    y(i+1,:) = y(i,:) + (9 * f(t(i+1), y(i+1,:)) + 19 * s(i,:) - 5 * s(i-1,:) + s(i-2,:)) * h24; % corrector
end;

Rewrite the code that counts the number of primes less than using an if-then statement.  Implement    your code where j goes from 2 to 15.

j:=1;
                               1
for i from 2 to 10
while ithprime(j)<2^i do
j:=j+1
end:
print(2^i,primes=j-1):
                        2048, primes = 9

I need to edit this code to satisfy a IF then Statement. can any one help me out?

\

regards "Geordi"

Hi everyone

I am currently trying to make my own simple package including a few procedures. So far I have been able to write some "code" that actually works when I open the document and hit "enter". I would, however, like to save the package so it can be accessed during any Maple session using the command "with". I have unsuccesfully tried to comprehend the Maple help pages regarding this question but I definitely don't want to mess things up.

This is what I have written:

mat := module ()
description "useful procedures for mathematics, physics and chemistry";
export AtomicWeight;
option package;

   AtomicWeight := proc (x) description "returns the average atomic mass of the naturally ocurring element";
   Units:-AddSystem(NewSystem, Units:-GetSystem(SI), u);
   return evalf(ScientificConstants:-Element(x, atomicweight, system = NewSystem, units))
   end proc

end module;

What should I do to save it correctly?

Thank in advance,
Mads


 

 

 

 

Following my previous question

http://www.mapleprimes.com/questions/200627-Lssolve-Midpoint

I wrote the following code

 

restart:
Phiavg:=0.06;
lambda:=0.05;
Ha:=0;
NBT:=0.5;
Nr:=500;
#N[bt]:=cc*NBT+(1-cc)*4; ## cc between 0 and 1
N[bt]:=cc*NBT+(1-cc^2)*0.75;


                              0.06
                              0.05
                               0
                              0.5
                              500
                                           2
                    0.5 cc + 0.75 - 0.75 cc
eq1:=diff(u(eta),eta,eta)+1/(mu(eta)/mu1[w])*(sigma-Nr*(phi(eta)-phi1[w])-(1-phi(eta))*T(eta)-Ha^2*u(eta))+((1/mu(eta)*(mu_phi*diff(phi(eta),eta)))*diff(u(eta),eta));
eq2:=diff(T(eta),eta)-1/(k(eta)/k1[w]);
eq3:=diff(phi(eta),eta)-phi(eta)/(N[bt]*(1-gama1*T(eta))^2)*diff(T(eta),eta);
 /  d   /  d         \\      1                                 
 |----- |----- u(eta)|| + ------- (mu1[w] (sigma - 500 phi(eta)
 \ deta \ deta       //   mu(eta)                              

    + 500 phi1[w] - (1 - phi(eta)) T(eta)))

             /  d           \ /  d         \
      mu_phi |----- phi(eta)| |----- u(eta)|
             \ deta         / \ deta       /
    + --------------------------------------
                     mu(eta)                
                    /  d         \   k1[w]
                    |----- T(eta)| - ------
                    \ deta       /   k(eta)
                                       /  d         \            
                              phi(eta) |----- T(eta)|            
/  d           \                       \ deta       /            
|----- phi(eta)| - ----------------------------------------------
\ deta         /   /                       2\                   2
                   \0.5 cc + 0.75 - 0.75 cc / (1 - gama1 T(eta))
mu:=unapply(mu1[bf]*(1+a[mu1]*phi(eta)+b[mu1]*phi(eta)^2),eta):
k:=unapply(k1[bf]*(1+a[k1]*phi(eta)+b[k1]*phi(eta)^2),eta):
rhop:=3880:
rhobf:=998.2:
cp:=773:
cbf:=4182:
rho:=unapply(  phi(eta)*rhop+(1-phi(eta))*rhobf ,eta):
c:=unapply(  (phi(eta)*rhop*cp+(1-phi(eta))*rhobf*cbf )/rho(eta) ,eta):
mu_phi:=mu1[bf]*(a[mu1]+2*b[mu1]*phi(eta)):
gama1:=0.00:
a[mu1]:=39.11:
b[mu1]:=533.9:
mu1[bf]:=9.93/10000:
a[k1]:=7.47:
b[k1]:=0:
k1[bf]:=0.597:
zet:=1:
phi1[w]:=phi0:
mu1[w]:=mu(0):
k1[w]:=k(0):

eq1:=subs(phi(0)=phi0,eq1);
eq2:=subs(phi(0)=phi0,eq2);
eq3:=subs(phi(0)=phi0,eq3);
/  d   /  d         \\   //                                    
|----- |----- u(eta)|| + \\0.0009930000000 + 0.03883623000 phi0
\ deta \ deta       //                                         

                      2\                                 
   + 0.5301627000 phi0 / (sigma - 500 phi(eta) + 500 phi0

                           \//               
   - (1 - phi(eta)) T(eta))/ \0.0009930000000

                                                   2\   
   + 0.03883623000 phi(eta) + 0.5301627000 phi(eta) / +

  /                                       /  d           \ /  d  
  |(0.03883623000 + 1.060325400 phi(eta)) |----- phi(eta)| |-----
  \                                       \ deta         / \ deta

         \\//                                        
   u(eta)|| \0.0009930000000 + 0.03883623000 phi(eta)
         //                                          

                          2\
   + 0.5301627000 phi(eta) /
           /  d         \     0.597 + 4.45959 phi0  
           |----- T(eta)| - ------------------------
           \ deta       /   0.597 + 4.45959 phi(eta)
                                        /  d         \
                            1. phi(eta) |----- T(eta)|
         /  d           \               \ deta       /
         |----- phi(eta)| - --------------------------
         \ deta         /                           2
                             0.5 cc + 0.75 - 0.75 cc  
Q:=proc(pp2,fi0) option remember; local res,F0,F1,F2,a,INT0,INT10,B;
print(pp2,fi0);
if not type([pp2,fi0],list(numeric)) then return 'procname(_passed)' end if;
res := dsolve(subs(sigma=pp2,phi0=fi0,{eq1=0,eq2=0,eq3=0,u(1)=-lambda*D(u)(1),u(0)=lambda*D(u)(0),phi(0)=phi0,T(0)=0}), numeric,output=listprocedure,initmesh=10, continuation=cc);
F0,F1,F2:=op(subs(res,[u(eta),phi(eta),T(eta)]));
INT0:=evalf(Int((abs(F0(eta)),eta=0..1)));
INT10:=evalf(Int(abs(F0(eta))*F1(eta),eta=0..1));
a[1]:=evalf(Int(F0(eta)*(F1(eta)*rhop+(1-F1(eta))*rhobf),eta=0..1));
#a[1]:=evalf(Int((F0(eta),eta=0..1)));
a[2]:=(INT10/INT0-Phiavg)/Phiavg; #relative
[a[1],a[2]]
end proc:
Q1:=proc(pp2,fi0) Q(_passed)[1] end proc;
Q2:=proc(pp2,fi0) Q(_passed)[2] end proc;
proc(pp2, fi0)  ...  end;
proc(pp2, fi0)  ...  end;
#Q(116,0.0041);
#tempe:=Optimization:-LSSolve([Q1,Q2],initialpoint=[130,0.01]);
#tempe:=Optimization:-LSSolve([Q1,Q2],initialpoint=[43.55,0.39]);
tempe:=Optimization:-LSSolve([Q1,Q2],initialpoint=[5.65,0.00036]);
#tempe:=Optimization:-LSSolve([Q1,Q2],initialpoint=[12,0.003]); # khoob ba 1
#tempe:=Optimization:-LSSolve([Q1,Q2],initialpoint=[5,0.01]);
                  HFloat(5.65), HFloat(3.6e-4)
           HFloat(5.650000070103341), HFloat(3.6e-4)
           HFloat(5.65), HFloat(3.600105456508193e-4)
     HFloat(29.63242379055208), HFloat(0.0205927592420527)
    HFloat(12.803902258015825), HFloat(0.006395385884750864)
    HFloat(12.803902403534572), HFloat(0.006395385884750864)
    HFloat(12.803902258015825), HFloat(0.00639539649402585)
   HFloat(12.804004931505949), HFloat(0.0063954867657199386)
    HFloat(12.804107604996073), HFloat(0.006395587646689013)
    HFloat(12.80400483062498), HFloat(0.006498160255844027)
    HFloat(12.803902157134855), HFloat(0.006498059374874952)
   HFloat(-1.0206939292143726), HFloat(-3.32764179807047e-4)
   HFloat(-1.0206939079125088), HFloat(-3.32764179807047e-4)
   HFloat(-1.0206939292143726), HFloat(-3.327536344433438e-4)
    HFloat(18.749500943683863), HFloat(0.01993840615828979)
    HFloat(3.9953780262640484), HFloat(0.00481041471606933)
     HFloat(6.166152606930136), HFloat(0.00703619658484674)
    HFloat(7.3193201827812295), HFloat(0.008218585352824569)
Error, (in Optimization:-LSSolve) complex value encountered
sigma:=tempe[2](1);
                          tempe[2](1)
phi0:=tempe[2](2);
                          tempe[2](2)
with(plots):

res2 := dsolve({eq1=0,eq2=0,eq3=0,u(1)=-lambda*D(u)(1),u(0)=lambda*D(u)(0),phi(0)=phi0,T(0)=0}, numeric,output=listprocedure,continuation=cc);
Error, (in dsolve/numeric/process_input) boundary conditions specified at too many points: {0, 1, 2}, can only solve two-point boundary value problems
G0,G1,G2:=op(subs(res2,[u(eta),phi(eta),T(eta)])):
ruu:=evalf((Int(abs(G0(eta))*(G1(eta)*rhop+(1-G1(eta))*rhobf ),eta=0..zet)))/(Phiavg*rhop+(1-Phiavg)*rhobf);
phb:=evalf((Int(abs(G0(eta))*G1(eta),eta=0..1))) / evalf((Int(abs(G0(eta)),eta=0..1))) ;
TTb:=evalf(Int(abs(G0(eta))*G2(eta)*(G1(eta)*rhop*cp+(1-G1(eta))*rhobf*cbf ),eta=0..1))/evalf(Int(abs(G0(eta))*(G1(eta)*rhop*cp+(1-G1(eta))*rhobf*cbf ),eta=0..1));
Error, invalid input: subs received res2, which is not valid for its 1st argument
                /  /1.                                        \
                | |                                           |
0.0008538922115 | |    |G0(eta)| (2881.8 G1(eta) + 998.2) deta|
                | |                                           |
                \/0.                                          /
                    /1.                       
                   |                          
                   |    |G0(eta)| G1(eta) deta
                   |                          
                  /0.                         
                  ----------------------------
                        /1.                   
                       |                      
                       |                      
                       |    |G0(eta)| deta    
                      /                       
                       0.                     
                                                              /Int(
                              1                               |     
------------------------------------------------------------- |     
  /1.                                                         |     
 |                                                            \     
 |              /             6                       6\            
 |    |G0(eta)| \-1.1752324 10  G1(eta) + 4.1744724 10 / deta       
/                                                                   
 0.                                                                 

                    /             6                       6\ , eta = 0. .. 1.)
  |G0(eta)| G2(eta) \-1.1752324 10  G1(eta) + 4.1744724 10 /                  

  \
  |
  |
  |
  /
#rhouu:=evalf((Int((G1(eta)*rhop+(1-G1(eta))*rhobf)*G0(eta),eta=0..1)));

odeplot(res2,[[eta,u(eta)/ruu],[eta,phi(eta)/phb],[eta,T(eta)/TTb]],0..1);
#odeplot(res2,[[eta,u(eta)],[eta,phi(eta)],[eta,T(eta)]],0..1);
Error, (in plots/odeplot) input is not a valid dsolve/numeric solution
rhou:=evalf((Int(abs(G0(eta))*(G1(eta)*rhop+(1-G1(eta))*rhobf ),eta=0..zet))):

Nub:=(1/G2(1))*(((1+a[k1]*abs(G1(0))+b[k1]*abs(G1(0))^2)/(1+a[k1]*Phiavg+b[k1]*Phiavg^2)));
                0.6905123602 (1 + 7.47 |G1(0)|)
                -------------------------------
                             G2(1)             
(rhs(res2(0.0000000000001)[3])-rhs(res2(0)[3]))/0.0000000000001;
Error, invalid input: rhs received res2(0.1e-12)[3], which is not valid for its 1st argument, expr
sigma;
                          tempe[2](1)
NBT;
                              0.5
>

 

the above code has been worked for NBT=0.6 and higher, whereas as NBT decreases, the code doesnt converge easily.

How can I fix this problem?

Thanks for your attention in advance

Amir

1 2 3 4 5 6 7 Last Page 1 of 10