Check if a given function is the solution of a PDE in Matlab

Question

Luca Donati 2023-5-26

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1974264-check-if-a-given-function-is-the-solution-of-a-pde-in-matlab

回答： Karan Singh 2023-8-29

在 MATLAB Online 中打开

I'm working on the following partial differential equation:

with

and

.

I have the solution which is:

with

and I'm trying to get the Green's function which should be:

with

and

. (

)

So to check that the Green function solves the differential equation in the forward variables y and s i wrote the following code:

clc
clear all
close all
a=0.3; mu=0.25; sigma=0.1; r0=0.1;
Mu = @(rt,s,t) mu+(rt-mu).*exp(-a.*(s-t))+((sigma^2)/(a^2)).*(exp(-a.*(s-t))-1);
Sigma = @(s,t) sqrt(sigma.^2./(2.*a) .* (1-exp(-2.*a.*(s-t))));
G= @(rs,s,rt,t) exp(-(rs-Mu(rt,s,t)).^2./(2.*Sigma(s,t).^2)) * 1./(sqrt(2*pi).*Sigma(s,t));
syms rs s rt t h;
DrsG=eval(['@(rs,s,rt,t)' char(diff(G,rs))]);
DrsrsG=eval(['@(rs,s,rt,t)' char(diff(DrsG,rs))]);
DsG=eval(['@(rs,s,rt,t)' char(diff(G,s))]);
Eq=@(rs,s,rt,t) DsG(rs,s,rt,t)+0.5*sigma^2*DrsrsG(rs,s,rt,t)+a*(mu-rs)*DrsG(rs,s,rt,t)-rs*G(rs,s,rt,t);
disp(Eq(0.019,1,0.01,0))

But it gives me -3.6936 (should be 0). So I thought maybe my green function is wrong and I wrote the following code to check if it gives me 0 for the correct solution (which I know it's correct):

F=@(h) (1-exp(-a.*(T-h)))./a;
B=@(h) (mu-sigma.^2./(2.*a.^2)).*(T-h-F(h)) + (F(h)).^2.*sigma.^2./(4.*a);
P=@(rt,h) exp(-rt.*F(h)-B(h));
DrtP=eval(['@(rt,h)' char(diff(P,rt))]);
DrtrtP=eval(['@(rt,h)' char(diff(DrtP,rt))]);
DhP=eval(['@(rt,h)' char(diff(P,h))]);
EqP=@(rt,h) DhP(rt,h)+0.5*sigma^2*DrtrtP(rt,h)+a*(mu-rt)*DrtP(rt,h)-rt*P(rt,h);
disp(EqP(0.019,1))

which gives me back -1.0408e-17. This is acceptably close to zero to consider the test correct right?

Is there any better way to check if the equation is varified? Also I got some values for a, mu and sigma from an article but maybe I could let them be syms and check everything in symbolic form?

4 个评论
显示 2更早的评论隐藏 2更早的评论

Torsten 2023-5-27

编辑：Torsten 2023-5-27

It's not safe that MATLAB at last will be able to satisfactory simplify the expressions you get so that it ends up with 0=0. But I would give it a try.

Luca Donati 2023-5-27

编辑：Luca Donati 2023-5-27

I tryed it but it isn't exact, I think it's because for

the Green's function isn't defined, so maybe there's a way to tell him to differentiate only for

?

Anyway using eval on the equation with the values I used on this code for the parameters and changing the variables, it's always very close to 0 so seems to be good.

请先登录，再进行评论。

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

Answer 1

Karan Singh 2023-8-29

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1974264-check-if-a-given-function-is-the-solution-of-a-pde-in-matlab#answer_1296026

在 MATLAB Online 中打开

Hi Luca, to verify if the equation is satisfied, you can check if the value is close to zero within a certain tolerance. In your case, the value “EqP(0.019,1)” is approximately “-1.0408e-17”, which is very close to zero. This suggests that the equation is indeed verified.

Keep in mind that numerical computations may introduce small errors due to finite precision arithmetic. Therefore, it is common to use a tolerance or threshold value to determine if a result is close enough to zero.

Regarding your question about using symbolic variables for “a”, “mu”, and “sigma”, it is indeed a good idea. By defining these variables as symbolic, you can perform symbolic computations and obtain exact solutions or expressions. This can help you in further analysis or manipulation of the equations.

To define these variables as symbolic, you can use the “syms ” function in MATLAB. Here's an example of how you can modify your code to use symbolic variables:

syms a mu sigma rs s rt t h;
Mu = mu + (rt - mu) * exp(-a * (s - t)) + ((sigma^2) / (a^2)) * (exp(-a * (s - t)) - 1);
Sigma = sqrt(sigma^2 / (2 * a) * (1 - exp(-2 * a * (s - t))));
G = exp(-(rs - Mu)^2 / (2 * Sigma^2)) * 1 / (sqrt(2 * pi) * Sigma);
DrsG = diff(G, rs);
DrsrsG = diff(DrsG, rs);
DsG = diff(G, s);
Eq = DsG + 0.5 * sigma^2 * DrsrsG + a * (mu - rs) * DrsG - rs * G;
disp(Eq);